### Hollibaugh and Rothenberg
### Agency Control through the Appointed Hierarchy: Presidential Politicization of Unilateral Appointees
### Main Code

rm(list=ls())

library(gtools)
library(lubridate)
library(plyr)
library(dplyr)
library(ggplot2)
library(nloptr)
library(scales)
library(lme4)
library(texreg)
library(MuMIn)
library(foreign)

### GLM Control function to assist in convergence
defaultControl <- list(algorithm="NLOPT_LN_BOBYQA",xtol_rel=1e-7,maxeval=1e6)
nloptwrap2 <- function(fn,par,lower,upper,control=list(),...) {
  for (n in names(defaultControl)) 
    if (is.null(control[[n]])) control[[n]] <- defaultControl[[n]]
    res <- nloptr(x0=par,eval_f=fn,lb=lower,ub=upper,opts=control,...)
    with(res,list(par=solution,
                  fval=objective,
                  feval=iterations,
                  conv=if (status>0) 0 else status,
                  message=message))
}

# mode utility function
Mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}


# asinh function to change scale on some plots
asinh_trans <- function(){
  trans_new(name = 'asinh', transform = function(x) asinh(x), 
            inverse = function(x) sinh(x))
}

# load in position file for Figure 1 
allposframe <- read.dta("HR-all-positions.dta")

# appropriate date format
allposframe$date <- as.Date(allposframe$date, format = "%Y-%m-%d", origin='1960-01-01')

# separate vline object for the graph
vline10 <- data.frame(type = "Noncareer SES Executives\n(As Percentage of Total SES Executives)",
                      vline = c(0.1, 0.25)) 

# figure 1
fig1 <- ggplot(allposframe, aes(x = date, y = positions, color = as.factor(agency != "All"), linetype = agency)) +
  geom_line() + 
  geom_line(data = subset(allposframe, agency == "All"), 
            aes(x = date, y = positions),
            linetype = 'solid',
            lwd = 1.25) + 
  scale_x_date(date_breaks = "2 years") + 
  geom_vline(xintercept = as.Date("1989-01-20"), linetype = "dashed") + 
  geom_vline(xintercept = as.Date("1993-01-20"), linetype = "dashed") + 
  geom_vline(xintercept = as.Date("2001-01-20"), linetype = "dashed") + 
  geom_vline(xintercept = as.Date("2009-01-20"), linetype = "dashed") + 
  geom_vline(xintercept = as.Date("2017-01-20"), linetype = "dashed") + 
  facet_wrap(~type, ncol = 1, scales = "free_y") + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "none") + 
  scale_color_grey(0.1, 0.9) + 
  scale_y_continuous(trans = 'log', breaks = c(0, 0.00001, 0.0001, 0.001, 0.01, 0.1)) + 
  xlab('\nDate') + 
  ylab('Proportion of Employees in Indicated Category\n') + 
  geom_hline(data = vline10, aes(yintercept = vline), linetype = "dotted")

ggsave(fig1, file = "Figure-1.pdf", height = 6, width = 8)
ggsave(fig1, file = "Figure-1.png", height = 6, width = 8)


rm(allposframe)

# load in main file
apptypes_merged <- read.dta("HR-appointments.dta")


# appropriate date format
apptypes_merged$date <- as.Date(apptypes_merged$date, format = "%Y-%m-%d", origin='1960-01-01')


# estimate a bunch of models!
schedc.glm.cong.meandur <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                     total - perm_excep_schedc - temp_excep_schedc) ~ 
                                 pres_switch + partyswitch + start_second_term + looming_election + 
                                 sameagency + z_pres_congdist * divided + z_polarization + 
                                 z_rollcalls + z_meanduration + prev_exec_pct + 
                                 factor(agency) + factor(pres_name),
                               data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                               family = quasibinomial)

schedc.glm.cong.simple.meandur <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                            total - perm_excep_schedc - temp_excep_schedc) ~ 
                                        z_pres_congdist * divided + z_meanduration + 
                                        factor(agency) + factor(pres_name),
                                      data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                                      family = quasibinomial)

schedc.glm.cong <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                             total - perm_excep_schedc - temp_excep_schedc) ~ 
                         pres_switch + partyswitch + start_second_term + looming_election + 
                         sameagency + z_pres_congdist * divided + z_polarization + 
                         z_rollcalls + prev_exec_pct + factor(agency) + factor(pres_name),
                       data = subset(apptypes_merged, agency != "All"), family = quasibinomial)

schedc.glm.cong.nodhs.nohhs <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                         total - perm_excep_schedc - temp_excep_schedc) ~ 
                                     pres_switch + partyswitch + start_second_term + looming_election + 
                                     sameagency + z_pres_congdist * divided + z_polarization + 
                                     z_rollcalls + prev_exec_pct + factor(agency) + factor(pres_name),
                                   data = subset(apptypes_merged, agency != "All" & agency != "HS" & agency != "HE"), 
                                   family = quasibinomial)

schedc.glm.cong.nodhs.nohhs.expertise <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                                   total - perm_excep_schedc - temp_excep_schedc) ~ 
                                               pres_switch + partyswitch + start_second_term + looming_election + 
                                               sameagency + z_pres_congdist * divided + z_polarization +
                                               z_rollcalls + prev_exec_pct + zmecompmedian + factor(agency) + factor(pres_name),
                                             data = subset(apptypes_merged, agency != "All" & agency != "HS" & agency != "HE"), 
                                             family = quasibinomial)

schedc.glm.fil.nodhs.nohhs <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                        total - perm_excep_schedc - temp_excep_schedc) ~ 
                                    pres_switch + partyswitch + start_second_term + looming_election + 
                                    sameagency + z_pres_fildist * divided + z_polarization + 
                                    z_rollcalls + prev_exec_pct + factor(agency) + factor(pres_name),
                                  data = subset(apptypes_merged, agency != "All" & agency != "HS" & agency != "HE"), 
                                  family = quasibinomial)

schedc.glm.fil.nodhs.nohhs.expertise <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                                  total - perm_excep_schedc - temp_excep_schedc) ~ 
                                              pres_switch + partyswitch + start_second_term + looming_election + 
                                              sameagency + z_pres_fildist * divided + z_polarization + 
                                              z_rollcalls + prev_exec_pct + zmecompmedian + factor(agency) + factor(pres_name),
                                            data = subset(apptypes_merged, agency != "All" & agency != "HS" & agency != "HE"), 
                                            family = quasibinomial)

schedc.glm.cong.nocsa <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                   total - perm_excep_schedc - temp_excep_schedc) ~ 
                               pres_switch + partyswitch + start_second_term + looming_election + 
                               sameagency + z_pres_congdist * divided + z_polarization + 
                               z_rollcalls + prev_exec_pct + factor(agency) + factor(pres_name),
                             data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS"))), 
                             family = quasibinomial)

schedc.glm.cong12 <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                               total12 - perm_excep_schedc - temp_excep_schedc) ~ 
                           pres_switch + partyswitch + start_second_term + looming_election + 
                           sameagency + z_pres_congdist * divided + z_polarization + 
                           z_rollcalls + prev_exec_pct + factor(agency) + factor(pres_name),
                         data = subset(apptypes_merged, agency != "All"), family = quasibinomial)

schedc.glm.cong12.expertise <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                         total12 - perm_excep_schedc - temp_excep_schedc) ~ 
                                     pres_switch + partyswitch + start_second_term + looming_election + 
                                     sameagency + z_pres_congdist * divided + z_polarization + 
                                     z_rollcalls + prev_exec_pct + zmecompmedian + factor(agency) + factor(pres_name),
                                   data = subset(apptypes_merged, agency != "All"), family = quasibinomial)

schedc.glm.cong.vacancies <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                             total - perm_excep_schedc - temp_excep_schedc) ~ 
                               pres_switch + partyswitch + start_second_term + looming_election + 
                               sameagency + z_pres_congdist * divided + z_polarization + 
                         z_rollcalls + prev_exec_pct + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                           factor(agency) + factor(pres_name),
                       data = subset(apptypes_merged, agency != "All"), family = quasibinomial)


schedc.glm.cong.vacancies.expertise <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                       total - perm_excep_schedc - temp_excep_schedc) ~ 
                                   pres_switch + partyswitch + start_second_term + looming_election + 
                                   sameagency + z_pres_congdist * divided + z_polarization + 
                                   z_rollcalls + prev_exec_pct + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                                   zmecompmedian + factor(agency) + factor(pres_name),
                                 data = subset(apptypes_merged, agency != "All"), family = quasibinomial)


schedc.glm.cong.expertise <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                       total - perm_excep_schedc - temp_excep_schedc) ~ 
                                   pres_switch + partyswitch + start_second_term + looming_election + 
                                   sameagency + z_pres_congdist * divided + z_polarization + 
                                   z_rollcalls + prev_exec_pct + zmecompmedian + 
                                   factor(agency) + factor(pres_name),
                                 data = subset(apptypes_merged, agency != "All"), family = quasibinomial)


schedc.glm.cong.vacancies.meandur <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                       total - perm_excep_schedc - temp_excep_schedc) ~ 
                                         pres_switch + partyswitch + start_second_term + looming_election + 
                                         sameagency + z_pres_congdist * divided + z_polarization + 
                                   z_rollcalls + z_meanduration + prev_exec_pct + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                                   factor(agency) + factor(pres_name),
                                 data = subset(apptypes_merged, agency != "All"), family = quasibinomial)

schedc.glm.cong.vacancies.expertise.meandur <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                               total - perm_excep_schedc - temp_excep_schedc) ~ 
                                           pres_switch + partyswitch + start_second_term + looming_election + 
                                           sameagency + z_pres_congdist * divided + z_polarization + 
                                           z_rollcalls + z_meanduration + prev_exec_pct + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                                           zmecompmedian + factor(agency) + factor(pres_name),
                                         data = subset(apptypes_merged, agency != "All"), family = quasibinomial)



schedc.glm.cong.expertise.meandur <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                       total - perm_excep_schedc - temp_excep_schedc) ~ 
                                         pres_switch + partyswitch + start_second_term + looming_election + 
                                         sameagency + z_pres_congdist * divided + z_polarization + 
                                   z_rollcalls + z_meanduration + prev_exec_pct + zmecompmedian + 
                                   factor(agency) + factor(pres_name),
                                 data = subset(apptypes_merged, agency != "All"), family = quasibinomial)


schedc.glm.cong.simple <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                             total - perm_excep_schedc - temp_excep_schedc) ~ 
                          z_pres_congdist * divided + factor(agency) + factor(pres_name),
                       data = subset(apptypes_merged, agency != "All"), family = quasibinomial)

schedc.glm.cong.simple.nodhs.nohhs <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                    total - perm_excep_schedc - temp_excep_schedc) ~ 
                                z_pres_congdist * divided + factor(agency) + factor(pres_name),
                              data = subset(apptypes_merged, agency != "All" & agency != "HS" & agency != "HE"), 
                              family = quasibinomial)

schedc.glm.fil.simple.nodhs.nohhs <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                                total - perm_excep_schedc - temp_excep_schedc) ~ 
                                            z_pres_fildist * divided + factor(agency) + factor(pres_name),
                                          data = subset(apptypes_merged, agency != "All" & agency != "HS" & agency != "HE"), 
                                          family = quasibinomial)


schedc.glm.cong.simple.nocsa <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                    total - perm_excep_schedc - temp_excep_schedc) ~ 
                                z_pres_congdist * divided + factor(agency) + factor(pres_name),
                              data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS"))), family = quasibinomial)


schedc.glm.cong12.simple <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                    total12 - perm_excep_schedc - temp_excep_schedc) ~ 
                                z_pres_congdist * divided + factor(agency) + factor(pres_name),
                              data = subset(apptypes_merged, agency != "All"), family = quasibinomial)



schedc.glmer.cong <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                 total - perm_excep_schedc - temp_excep_schedc) ~ 
                             pres_switch + partyswitch + start_second_term + looming_election + 
                             sameagency + z_pres_congdist * divided + z_polarization +
                             z_rollcalls + prev_exec_pct + central_staff_agency + (1|agency) + (1|pres_name),
                           data = subset(apptypes_merged, agency != "All"), family = binomial,
                           control=glmerControl(optimizer = "nloptwrap2",
                                                optCtrl=list(maxfun=2e4)))


schedc.glmer.cong.nodhs.nohhs <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                 total - perm_excep_schedc - temp_excep_schedc) ~ 
                                   pres_switch + partyswitch + start_second_term + looming_election + 
                                   sameagency + z_pres_congdist * divided + z_polarization +
                             z_rollcalls + prev_exec_pct + central_staff_agency + (1|agency) + (1|pres_name),
                           data = subset(apptypes_merged, agency != "All" & agency != "HS" & agency != "HE"), 
                           family = binomial,
                           control=glmerControl(optimizer = "nloptwrap2",
                                                optCtrl=list(maxfun=2e4)))


schedc.glmer.cong.nodhs.nohhs.expertise <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                             total - perm_excep_schedc - temp_excep_schedc) ~ 
                                         pres_switch + partyswitch + start_second_term + looming_election + 
                                         sameagency + z_pres_congdist * divided + z_polarization +
                                         z_rollcalls + prev_exec_pct + central_staff_agency + zmecompmedian + (1|agency) + (1|pres_name),
                                       data = subset(apptypes_merged, agency != "All" & agency != "HS" & agency != "HE"), 
                                       family = binomial,
                                       control=glmerControl(optimizer = "nloptwrap2",
                                                            optCtrl=list(maxfun=2e4)))

schedc.glmer.fil.nodhs.nohhs <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                             total - perm_excep_schedc - temp_excep_schedc) ~ 
                                        pres_switch + partyswitch + start_second_term + looming_election + 
                                        sameagency + z_pres_fildist * divided + z_polarization +
                                         z_rollcalls + prev_exec_pct + central_staff_agency + (1|agency) + (1|pres_name),
                                       data = subset(apptypes_merged, agency != "All" & agency != "HS" & agency != "HE"), 
                                       family = binomial,
                                       control=glmerControl(optimizer = "nloptwrap2",
                                                            optCtrl=list(maxfun=2e4)))


schedc.glmer.fil.nodhs.nohhs.expertise <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                            total - perm_excep_schedc - temp_excep_schedc) ~ 
                                        pres_switch + partyswitch + start_second_term + looming_election + 
                                        sameagency + z_pres_fildist * divided + z_polarization +
                                        z_rollcalls + prev_exec_pct + central_staff_agency + zmecompmedian + (1|agency) + (1|pres_name),
                                      data = subset(apptypes_merged, agency != "All" & agency != "HS" & agency != "HE"), 
                                      family = binomial,
                                      control=glmerControl(optimizer = "nloptwrap2",
                                                           optCtrl=list(maxfun=2e4)))


schedc.glmer.cong.nocsa <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                 total - perm_excep_schedc - temp_excep_schedc) ~ 
                                   pres_switch + partyswitch + start_second_term + looming_election + 
                                   sameagency + z_pres_congdist * divided + z_polarization +
                             z_rollcalls + prev_exec_pct + (1|agency) + (1|pres_name),
                           data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS"))), family = binomial,
                           control=glmerControl(optimizer = "nloptwrap2",
                                                optCtrl=list(maxfun=2e4)))


schedc.glmer.cong.nocsa.expertise <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                       total - perm_excep_schedc - temp_excep_schedc) ~ 
                                   pres_switch + partyswitch + start_second_term + looming_election + 
                                   sameagency + z_pres_congdist * divided + z_polarization +
                                   z_rollcalls + prev_exec_pct + zmecompmedian + (1|agency) + (1|pres_name),
                                 data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS"))), family = binomial,
                                 control=glmerControl(optimizer = "nloptwrap2",
                                                      optCtrl=list(maxfun=2e4)))


schedc.glmer.cong12 <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                 total12 - perm_excep_schedc - temp_excep_schedc) ~ 
                               pres_switch + partyswitch + start_second_term + looming_election + 
                               sameagency + z_pres_congdist * divided + z_polarization +
                             z_rollcalls + prev_exec_pct + central_staff_agency + (1|agency) + (1|pres_name),
                           data = subset(apptypes_merged, agency != "All"), family = binomial,
                           control=glmerControl(optimizer = "nloptwrap2",
                                                optCtrl=list(maxfun=2e4)))


schedc.glmer.cong12.expertise <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                   total12 - perm_excep_schedc - temp_excep_schedc) ~ 
                               pres_switch + partyswitch + start_second_term + looming_election + 
                               sameagency + z_pres_congdist * divided + z_polarization +
                               z_rollcalls + prev_exec_pct + central_staff_agency + zmecompmedian + (1|agency) + (1|pres_name),
                             data = subset(apptypes_merged, agency != "All"), family = binomial,
                             control=glmerControl(optimizer = "nloptwrap2",
                                                  optCtrl=list(maxfun=2e4)))


schedc.glmer.cong.vacancies <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                 total - perm_excep_schedc - temp_excep_schedc) ~ 
                                   pres_switch + partyswitch + start_second_term + looming_election + 
                                   sameagency + z_pres_congdist * divided + z_polarization +
                             z_rollcalls + prev_exec_pct +  I(prev_empty/(prev_empty + prev_interim + prev_permanent))  + 
                               (1|agency) + (1|pres_name),
                           data = subset(apptypes_merged, agency != "All"), family = binomial,
                           control=glmerControl(optimizer = "nloptwrap2",
                                                optCtrl=list(maxfun=2e4)))

schedc.glmer.cong.vacancies.expertise <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                           total - perm_excep_schedc - temp_excep_schedc) ~ 
                                       pres_switch + partyswitch + start_second_term + looming_election + 
                                       sameagency + z_pres_congdist * divided + z_polarization +
                                       z_rollcalls + prev_exec_pct +  I(prev_empty/(prev_empty + prev_interim + prev_permanent))  + 
                                       zmecompmedian + (1|agency) + (1|pres_name),
                                     data = subset(apptypes_merged, agency != "All"), family = binomial,
                                     control=glmerControl(optimizer = "nloptwrap2",
                                                          optCtrl=list(maxfun=2e4)))


schedc.glmer.cong.expertise <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                           total - perm_excep_schedc - temp_excep_schedc) ~ 
                                       pres_switch + partyswitch + start_second_term + looming_election + 
                                       sameagency + z_pres_congdist * divided + z_polarization +
                                       z_rollcalls + prev_exec_pct + central_staff_agency + zmecompmedian + 
                                       (1|agency) + (1|pres_name),
                                     data = subset(apptypes_merged, agency != "All"), family = binomial,
                                     control=glmerControl(optimizer = "nloptwrap2",
                                                          optCtrl=list(maxfun=2e4)))



schedc.glmer.cong.vacancies.meandur <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                           total - perm_excep_schedc - temp_excep_schedc) ~ 
                                             pres_switch + partyswitch + start_second_term + looming_election + 
                                             sameagency + z_pres_congdist * divided + z_polarization +
                                       z_rollcalls + z_meanduration + prev_exec_pct +  I(prev_empty/(prev_empty + prev_interim + prev_permanent))  + 
                                       (1|agency) + (1|pres_name),
                                     data = subset(apptypes_merged, agency != "All"), family = binomial,
                                     control=glmerControl(optimizer = "nloptwrap2",
                                                          optCtrl=list(maxfun=2e4)))

schedc.glmer.cong.vacancies.expertise.meandur <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                                   total - perm_excep_schedc - temp_excep_schedc) ~ 
                                               pres_switch + partyswitch + start_second_term + looming_election + 
                                               sameagency + z_pres_congdist * divided + z_polarization +
                                               z_rollcalls + z_meanduration + prev_exec_pct +  I(prev_empty/(prev_empty + prev_interim + prev_permanent))  + 
                                               zmecompmedian + (1|agency) + (1|pres_name),
                                             data = subset(apptypes_merged, agency != "All"), family = binomial,
                                             control=glmerControl(optimizer = "nloptwrap2",
                                                                  optCtrl=list(maxfun=2e4)))


schedc.glmer.cong.expertise.meandur <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                           total - perm_excep_schedc - temp_excep_schedc) ~ 
                                             pres_switch + partyswitch + start_second_term + looming_election + 
                                             sameagency + z_pres_congdist * divided + z_polarization +
                                       z_rollcalls + z_meanduration + prev_exec_pct + central_staff_agency + zmecompmedian + 
                                       (1|agency) + (1|pres_name),
                                     data = subset(apptypes_merged, agency != "All"), family = binomial,
                                     control=glmerControl(optimizer = "nloptwrap2",
                                                          optCtrl=list(maxfun=2e4)))



schedc.glmer.cong.def <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                 total - perm_excep_schedc - temp_excep_schedc) ~ 
                                 pres_switch + partyswitch + start_second_term + looming_election + 
                                 sameagency + z_pres_congdist * divided + z_polarization +
                             z_rollcalls + prev_exec_pct + (1|agency) + (1|pres_name),
                           data = subset(apptypes_merged, agency %in% c("DD", "AR", "AF", "NV")), family = binomial,
                           control=glmerControl(optimizer = "nloptwrap2",
                                                optCtrl=list(maxfun=2e4)))


schedc.glmer.cong.def.expertise <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                     total - perm_excep_schedc - temp_excep_schedc) ~ 
                                 pres_switch + partyswitch + start_second_term + looming_election + 
                                 sameagency + z_pres_congdist * divided + z_polarization +
                                 z_rollcalls + prev_exec_pct + zmecompmedian + (1|pres_name),
                               data = subset(apptypes_merged, agency %in% c("DD", "AR", "AF", "NV")), family = binomial,
                               control=glmerControl(optimizer = "nloptwrap2",
                                                    optCtrl=list(maxfun=2e4)))


schedc.glm.cong.def <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                     total - perm_excep_schedc - temp_excep_schedc) ~ 
                             pres_switch + partyswitch + start_second_term + looming_election + 
                             sameagency + z_pres_congdist * divided + z_polarization +
                                 z_rollcalls + prev_exec_pct + factor(agency) + factor(pres_name),
                               data = subset(apptypes_merged, agency %in% c("DD", "AR", "AF", "NV")), family = quasibinomial)

schedc.glmer.cong.civ <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                     total - perm_excep_schedc - temp_excep_schedc) ~ 
                                 pres_switch + partyswitch + start_second_term + looming_election + 
                                 sameagency + z_pres_congdist * divided + z_polarization +
                                 z_rollcalls + prev_exec_pct + central_staff_agency + (1|agency) + (1|pres_name),
                               data = subset(apptypes_merged, !(agency %in% c("All", "DD", "AR", "AF", "NV"))), family = binomial,
                               control=glmerControl(optimizer = "nloptwrap2",
                                                    optCtrl=list(maxfun=2e4)))

schedc.glmer.cong.civ.expertise <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                     total - perm_excep_schedc - temp_excep_schedc) ~ 
                                 pres_switch + partyswitch + start_second_term + looming_election + 
                                 sameagency + z_pres_congdist * divided + z_polarization +
                                 z_rollcalls + prev_exec_pct + central_staff_agency + zmecompmedian + (1|agency) + (1|pres_name),
                               data = subset(apptypes_merged, !(agency %in% c("All", "DD", "AR", "AF", "NV"))), family = binomial,
                               control=glmerControl(optimizer = "nloptwrap2",
                                                    optCtrl=list(maxfun=2e4)))


schedc.glmer.cong.meandur <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                 total - perm_excep_schedc - temp_excep_schedc) ~ 
                                   pres_switch + partyswitch + start_second_term + looming_election + 
                                   sameagency + z_pres_congdist * divided + z_polarization +
                             z_rollcalls + z_meanduration + prev_exec_pct + (1|agency) + (1|pres_name),
                           data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                           family = binomial,
                           control=glmerControl(optimizer = "nloptwrap2",
                                                optCtrl=list(maxfun=2e4)))

schedc.glmer.cong.simple <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                 total - perm_excep_schedc - temp_excep_schedc) ~ 
                             z_pres_congdist * divided + (1|agency) + (1|pres_name),
                           data = subset(apptypes_merged, agency != "All"), family = binomial,
                           control=glmerControl(optimizer = "nloptwrap2",
                                                optCtrl=list(maxfun=2e4)))


schedc.glmer.cong.simple.nodhs.nohhs <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                        total - perm_excep_schedc - temp_excep_schedc) ~ 
                                    z_pres_congdist * divided + (1|agency) + (1|pres_name),
                                  data = subset(apptypes_merged, agency != "All" & agency != "HS" & agency != "HE"), 
                                  family = binomial,
                                  control=glmerControl(optimizer = "nloptwrap2",
                                                       optCtrl=list(maxfun=2e4)))


schedc.glmer.fil.simple.nodhs.nohhs <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                                    total - perm_excep_schedc - temp_excep_schedc) ~ 
                                                z_pres_fildist * divided + (1|agency) + (1|pres_name),
                                              data = subset(apptypes_merged, agency != "All" & agency != "HS" & agency != "HE"), 
                                              family = binomial,
                                              control=glmerControl(optimizer = "nloptwrap2",
                                                                   optCtrl=list(maxfun=2e4)))



schedc.glmer.cong.simple.nocsa <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                        total - perm_excep_schedc - temp_excep_schedc) ~ 
                                    z_pres_congdist * divided + (1|agency) + (1|pres_name),
                                  data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS"))), family = binomial,
                                  control=glmerControl(optimizer = "nloptwrap2",
                                                       optCtrl=list(maxfun=2e4)))



schedc.glmer.cong12.simple <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                        total12 - perm_excep_schedc - temp_excep_schedc) ~ 
                                    z_pres_congdist * divided + (1|agency) + (1|pres_name),
                                  data = subset(apptypes_merged, agency != "All"), family = binomial,
                                  control=glmerControl(optimizer = "nloptwrap2",
                                                       optCtrl=list(maxfun=2e4)))


schedc.glmer.cong.simple.def <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                        total - perm_excep_schedc - temp_excep_schedc) ~ 
                                    z_pres_congdist * divided + (1|agency) + (1|pres_name),
                                  data = subset(apptypes_merged, agency %in% c("DD", "AF", "AR", "NV")), family = binomial,
                                  control=glmerControl(optimizer = "nloptwrap2",
                                                       optCtrl=list(maxfun=2e4)))

schedc.glm.cong.simple.def <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                            total - perm_excep_schedc - temp_excep_schedc) ~ 
                                        z_pres_congdist * divided + factor(agency) + factor(pres_name),
                                      data = subset(apptypes_merged, agency %in% c("DD", "AF", "AR", "NV")), family = quasibinomial)


schedc.glmer.cong.simple.civ <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                            total - perm_excep_schedc - temp_excep_schedc) ~ 
                                        z_pres_congdist * divided + (1|agency) + (1|pres_name),
                                      data = subset(apptypes_merged, !(agency %in% c("All", "DD", "AF", "AR", "NV"))), family = binomial,
                                      control=glmerControl(optimizer = "nloptwrap2",
                                                           optCtrl=list(maxfun=2e4)))

schedc.glmer.cong.only.meandur <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                                total - perm_excep_schedc - temp_excep_schedc) ~ 
                                            z_meanduration + (1|agency) + (1|pres_name),
                                        data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                                        family = binomial,
                                          control=glmerControl(optimizer = "nloptwrap2",
                                                               optCtrl=list(maxfun=2e4)))


schedc.glm.cong.only.meandur <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                              total - perm_excep_schedc - temp_excep_schedc) ~ 
                                          z_meanduration + factor(agency) + factor(pres_name),
                                        data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                                        family = quasibinomial)


schedc.glmer.cong.simple.meandur <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                        total - perm_excep_schedc - temp_excep_schedc) ~ 
                                    z_pres_congdist * divided + z_meanduration + (1|agency) + (1|pres_name),
                                    data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                                    family = binomial,
                                  control=glmerControl(optimizer = "nloptwrap2",
                                                       optCtrl=list(maxfun=2e4)))


schedc.glm.cong.cj.meandur <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                total - perm_excep_schedc - temp_excep_schedc) ~ 
                            partyswitch + start_second_term + looming_election + 
                            z_pres_agencydist + z_pres_congdist * divided + z_polarization +
                            z_rollcalls + z_meanduration + prev_exec_pct + factor(agency) + factor(pres_name),
                            data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                            family = quasibinomial)

schedc.glm.cong.cj.expertise.meandur <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                        total - perm_excep_schedc - temp_excep_schedc) ~ 
                                    partyswitch + start_second_term + looming_election + 
                                    z_pres_agencydist + z_pres_congdist * divided + z_polarization +
                                    z_rollcalls + z_meanduration + prev_exec_pct + zmecompmedian + factor(agency) + factor(pres_name),
                                  data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                                  family = quasibinomial)


schedc.glm.cong.cj.expertise.meandur <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                                  total - perm_excep_schedc - temp_excep_schedc) ~ 
                                              partyswitch + start_second_term + looming_election + 
                                              z_pres_agencydist + z_pres_congdist * divided + z_polarization +
                                              z_rollcalls + z_meanduration + prev_exec_pct + zmecompmedian + factor(agency) + factor(pres_name),
                                            data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                                            family = quasibinomial)


schedc.glm.cong.cj <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                             total - perm_excep_schedc - temp_excep_schedc) ~ 
                            partyswitch + start_second_term + looming_election + 
                           z_pres_agencydist + z_pres_congdist * divided + z_polarization +
                           z_rollcalls + prev_exec_pct + factor(agency) + factor(pres_name),
                       data = subset(apptypes_merged, agency != "All"), family = quasibinomial)

schedc.glm.cong.cj.expertise <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                total - perm_excep_schedc - temp_excep_schedc) ~ 
                            partyswitch + start_second_term + looming_election + 
                            z_pres_agencydist + z_pres_congdist * divided + z_polarization +
                            z_rollcalls + prev_exec_pct + zmecompmedian + factor(agency) + factor(pres_name),
                          data = subset(apptypes_merged, agency != "All"), family = quasibinomial)


schedc.glm.fil.cj <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                total - perm_excep_schedc - temp_excep_schedc) ~ 
                            partyswitch + start_second_term + looming_election + 
                            z_pres_agencydist + z_pres_fildist * divided + z_polarization +
                            z_rollcalls + prev_exec_pct + factor(agency) + factor(pres_name),
                          data = subset(apptypes_merged, agency != "All"), family = quasibinomial)


schedc.glmer.cong.cj.meandur <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                    total - perm_excep_schedc - temp_excep_schedc) ~ 
                                partyswitch + start_second_term + looming_election + 
                                z_pres_agencydist + z_pres_congdist * divided + z_polarization +
                                z_rollcalls + z_meanduration + prev_exec_pct + (1|agency) + (1|pres_name),
                                data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                                family = binomial,
                              control=glmerControl(optimizer = "nloptwrap2",
                                                   optCtrl=list(maxfun=2e4)))


schedc.glmer.cong.cj.expertise.meandur <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                            total - perm_excep_schedc - temp_excep_schedc) ~ 
                                        partyswitch + start_second_term + looming_election + 
                                        z_pres_agencydist + z_pres_congdist * divided + z_polarization +
                                        z_rollcalls + z_meanduration + prev_exec_pct + zmecompmedian + (1|agency) + (1|pres_name),
                                      data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                                      family = binomial,
                                      control=glmerControl(optimizer = "nloptwrap2",
                                                           optCtrl=list(maxfun=2e4)))

schedc.glmer.cong.cj <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                total - perm_excep_schedc - temp_excep_schedc) ~ 
                                partyswitch + start_second_term + looming_election + 
                               z_pres_agencydist + z_pres_congdist * divided + z_polarization +
                             z_rollcalls + central_staff_agency + prev_exec_pct + (1|agency) + (1|pres_name),
                          data = subset(apptypes_merged, agency != "All"), family = binomial,
                          control=glmerControl(optimizer = "nloptwrap2",
                                               optCtrl=list(maxfun=2e4)))

schedc.glmer.cong.cj.expertise <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                    total - perm_excep_schedc - temp_excep_schedc) ~ 
                                partyswitch + start_second_term + looming_election + 
                                z_pres_agencydist + z_pres_congdist * divided + z_polarization +
                                z_rollcalls + central_staff_agency + prev_exec_pct + zmecompmedian + (1|agency) + (1|pres_name),
                              data = subset(apptypes_merged, agency != "All"), family = binomial,
                              control=glmerControl(optimizer = "nloptwrap2",
                                                   optCtrl=list(maxfun=2e4)))


schedc.glm.fil.simple <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                            total - perm_excep_schedc - temp_excep_schedc) ~ 
                        z_pres_fildist * divided + factor(agency) + factor(pres_name),
                      data = subset(apptypes_merged, agency != "All"), family = quasibinomial)

schedc.glm.fil.simple.meandur <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                            total - perm_excep_schedc - temp_excep_schedc) ~ 
                                       z_pres_fildist * divided + z_meanduration + factor(agency) + factor(pres_name),
                                      data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                                      family = quasibinomial)


schedc.glm.fil.simple.nocsa <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                   total - perm_excep_schedc - temp_excep_schedc) ~ 
                               z_pres_fildist * divided + factor(agency) + factor(pres_name),
                             data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS"))), family = quasibinomial)


schedc.glm.fil12.simple <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                   total12 - perm_excep_schedc - temp_excep_schedc) ~ 
                               z_pres_fildist * divided + factor(agency) + factor(pres_name),
                             data = subset(apptypes_merged, agency != "All"), family = quasibinomial)


schedc.glm.fil.meandur <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                            total - perm_excep_schedc - temp_excep_schedc) ~ 
                              pres_switch + partyswitch + start_second_term + looming_election + 
                              sameagency + z_pres_fildist * divided + z_polarization +
                        z_rollcalls + z_meanduration + prev_exec_pct + factor(agency) + factor(pres_name),
                        data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                        family = quasibinomial)

schedc.glmer.fil.only.meandur <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                               total - perm_excep_schedc - temp_excep_schedc) ~ 
                                           z_meanduration + (1|agency) + (1|pres_name),
                                       data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                                       family = binomial,
                                         control=glmerControl(optimizer = "nloptwrap2",
                                                              optCtrl=list(maxfun=2e4)))

schedc.glm.fil.only.meandur <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                             total - perm_excep_schedc - temp_excep_schedc) ~ 
                                         z_meanduration + factor(agency) + factor(pres_name),
                                       data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                                       family = quasibinomial)



schedc.glmer.fil.simple.meandur <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                                total - perm_excep_schedc - temp_excep_schedc) ~ 
                                            z_pres_fildist * divided + z_meanduration + (1|agency) + (1|pres_name),
                                         data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                                         family = binomial,
                                          control=glmerControl(optimizer = "nloptwrap2",
                                                               optCtrl=list(maxfun=2e4)))

schedc.glm.fil <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                total - perm_excep_schedc - temp_excep_schedc) ~ 
                        pres_switch + partyswitch + start_second_term + looming_election + 
                        sameagency + z_pres_fildist * divided + z_polarization +
                            z_rollcalls + prev_exec_pct  + factor(agency) + factor(pres_name),
                          data = subset(apptypes_merged, agency != "All"), family = quasibinomial)


schedc.glm.fil.nocsa <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                            total - perm_excep_schedc - temp_excep_schedc) ~ 
                              pres_switch + partyswitch + start_second_term + looming_election + 
                              sameagency + z_pres_fildist * divided + z_polarization +
                        z_rollcalls + prev_exec_pct  + factor(agency) + factor(pres_name),
                      data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS"))), family = quasibinomial)


schedc.glm.fil12 <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                            total12 - perm_excep_schedc - temp_excep_schedc) ~ 
                          pres_switch + partyswitch + start_second_term + looming_election + 
                          sameagency + z_pres_fildist * divided + z_polarization +
                        z_rollcalls + prev_exec_pct  + factor(agency) + factor(pres_name),
                      data = subset(apptypes_merged, agency != "All"), family = quasibinomial)


schedc.glm.fil12.expertise <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                              total12 - perm_excep_schedc - temp_excep_schedc) ~ 
                          pres_switch + partyswitch + start_second_term + looming_election + 
                          sameagency + z_pres_fildist * divided + z_polarization +
                          z_rollcalls + prev_exec_pct  + zmecompmedian + factor(agency) + factor(pres_name),
                        data = subset(apptypes_merged, agency != "All"), family = quasibinomial)

schedc.glm.fil.vacancies <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                            total - perm_excep_schedc - temp_excep_schedc) ~ 
                              pres_switch + partyswitch + start_second_term + looming_election + 
                              sameagency + z_pres_fildist * divided + z_polarization +
                        z_rollcalls + prev_exec_pct  + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                          factor(agency) + factor(pres_name),
                      data = subset(apptypes_merged, agency != "All"), family = quasibinomial)

schedc.glm.fil.vacancies.expertise <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                      total - perm_excep_schedc - temp_excep_schedc) ~ 
                                  pres_switch + partyswitch + start_second_term + looming_election + 
                                  sameagency + z_pres_fildist * divided + z_polarization +
                                  z_rollcalls + prev_exec_pct  + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                                  zmecompmedian + factor(agency) + factor(pres_name),
                                data = subset(apptypes_merged, agency != "All"), family = quasibinomial)


schedc.glm.fil.expertise <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                      total - perm_excep_schedc - temp_excep_schedc) ~ 
                                  pres_switch + partyswitch + start_second_term + looming_election + 
                                  sameagency + z_pres_fildist * divided + z_polarization +
                                  z_rollcalls + prev_exec_pct  + zmecompmedian + 
                                  factor(agency) + factor(pres_name),
                                data = subset(apptypes_merged, agency != "All"), family = quasibinomial)

schedc.glm.fil.vacancies.meandur <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                      total - perm_excep_schedc - temp_excep_schedc) ~ 
                                        pres_switch + partyswitch + start_second_term + looming_election + 
                                        sameagency + z_pres_fildist * divided + z_polarization +
                                  z_rollcalls + z_meanduration + prev_exec_pct  + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                                  factor(agency) + factor(pres_name),
                                data = subset(apptypes_merged, agency != "All"), family = quasibinomial)

schedc.glm.fil.vacancies.expertise.meandur <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                              total - perm_excep_schedc - temp_excep_schedc) ~ 
                                          pres_switch + partyswitch + start_second_term + looming_election + 
                                          sameagency + z_pres_fildist * divided + z_polarization +
                                          z_rollcalls + z_meanduration + prev_exec_pct  + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                                          zmecompmedian + factor(agency) + factor(pres_name),
                                        data = subset(apptypes_merged, agency != "All"), family = quasibinomial)


schedc.glm.fil.expertise.meandur <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                      total - perm_excep_schedc - temp_excep_schedc) ~ 
                                        pres_switch + partyswitch + start_second_term + looming_election + 
                                        sameagency + z_pres_fildist * divided + z_polarization +
                                  z_rollcalls + z_meanduration + prev_exec_pct + zmecompmedian + 
                                  factor(agency) + factor(pres_name),
                                data = subset(apptypes_merged, agency != "All"), family = quasibinomial)


schedc.glmer.fil.simple <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                total - perm_excep_schedc - temp_excep_schedc) ~ 
                            z_pres_fildist * divided + 
                              (1|agency) + (1|pres_name),
                          data = subset(apptypes_merged, agency != "All"), family = binomial,
                          control=glmerControl(optimizer = "nloptwrap2",
                                               optCtrl=list(maxfun=2e4)))

schedc.glmer.fil.simple.nocsa <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                       total - perm_excep_schedc - temp_excep_schedc) ~ 
                                   z_pres_fildist * divided + 
                                   (1|agency) + (1|pres_name),
                                 data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS"))), family = binomial,
                                 control=glmerControl(optimizer = "nloptwrap2",
                                                      optCtrl=list(maxfun=2e4)))


schedc.glmer.fil12.simple <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                       total12 - perm_excep_schedc - temp_excep_schedc) ~ 
                                   z_pres_fildist * divided + 
                                   (1|agency) + (1|pres_name),
                                 data = subset(apptypes_merged, agency != "All"), family = binomial,
                                 control=glmerControl(optimizer = "nloptwrap2",
                                                      optCtrl=list(maxfun=2e4)))


schedc.glmer.fil.simple.def <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                       total - perm_excep_schedc - temp_excep_schedc) ~ 
                                   z_pres_fildist * divided + (1|agency) + (1|pres_name),
                                 data = subset(apptypes_merged, agency %in% c("DD", "AF", "AR", "NV")), family = binomial,
                                 control=glmerControl(optimizer = "nloptwrap2",
                                                      optCtrl=list(maxfun=2e4)))

schedc.glm.fil.simple.def <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                           total - perm_excep_schedc - temp_excep_schedc) ~ 
                                       z_pres_fildist * divided + factor(agency) + factor(pres_name),
                                     data = subset(apptypes_merged, agency %in% c("DD", "AF", "AR", "NV")), family = quasibinomial)


schedc.glmer.fil.simple.civ <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                           total - perm_excep_schedc - temp_excep_schedc) ~ 
                                       z_pres_fildist * divided + (1|agency) + (1|pres_name),
                                     data = subset(apptypes_merged, !(agency %in% c("All", "DD", "AF", "AR", "NV"))), family = binomial,
                                     control=glmerControl(optimizer = "nloptwrap2",
                                                          optCtrl=list(maxfun=2e4)))

schedc.glmer.fil <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                total - perm_excep_schedc - temp_excep_schedc) ~ 
                            pres_switch + partyswitch + start_second_term + looming_election + 
                            sameagency + z_pres_fildist * divided + z_polarization +
                            z_rollcalls +  prev_exec_pct + central_staff_agency + (1|agency) + (1|pres_name),
                          data = subset(apptypes_merged, agency != "All"), family = binomial,
                          control=glmerControl(optimizer = "nloptwrap2",
                                               optCtrl=list(maxfun=2e4)))


schedc.glmer.fil.nocsa <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                total - perm_excep_schedc - temp_excep_schedc) ~ 
                                  pres_switch + partyswitch + start_second_term + looming_election + 
                                  sameagency + z_pres_fildist * divided + z_polarization +
                            z_rollcalls +  prev_exec_pct + central_staff_agency + (1|agency) + (1|pres_name),
                          data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS"))), family = binomial,
                          control=glmerControl(optimizer = "nloptwrap2",
                                               optCtrl=list(maxfun=2e4)))


schedc.glmer.fil.nocsa.expertise <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                      total - perm_excep_schedc - temp_excep_schedc) ~ 
                                  pres_switch + partyswitch + start_second_term + looming_election + 
                                  sameagency + z_pres_fildist * divided + z_polarization +
                                  z_rollcalls +  prev_exec_pct + zmecompmedian + (1|agency) + (1|pres_name),
                                data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS"))), family = binomial,
                                control=glmerControl(optimizer = "nloptwrap2",
                                                     optCtrl=list(maxfun=2e4)))

schedc.glmer.fil12 <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                total12 - perm_excep_schedc - temp_excep_schedc) ~ 
                              pres_switch + partyswitch + start_second_term + looming_election + 
                              sameagency + z_pres_fildist * divided + z_polarization +
                            z_rollcalls +  prev_exec_pct + central_staff_agency + (1|agency) + (1|pres_name),
                          data = subset(apptypes_merged, agency != "All"), family = binomial,
                          control=glmerControl(optimizer = "nloptwrap2",
                                               optCtrl=list(maxfun=2e4)))

schedc.glmer.fil12.expertise <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                  total12 - perm_excep_schedc - temp_excep_schedc) ~ 
                              pres_switch + partyswitch + start_second_term + looming_election + 
                              sameagency + z_pres_fildist * divided + z_polarization +
                              z_rollcalls +  prev_exec_pct + central_staff_agency + zmecompmedian + (1|agency) + (1|pres_name),
                            data = subset(apptypes_merged, agency != "All"), family = binomial,
                            control=glmerControl(optimizer = "nloptwrap2",
                                                 optCtrl=list(maxfun=2e4)))



schedc.glmer.fil.vacancies <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                total - perm_excep_schedc - temp_excep_schedc) ~ 
                                  pres_switch + partyswitch + start_second_term + looming_election + 
                                  sameagency + z_pres_fildist * divided + z_polarization +
                            z_rollcalls +  prev_exec_pct + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                            (1|agency) + (1|pres_name),
                          data = subset(apptypes_merged, agency != "All"), family = binomial,
                          control=glmerControl(optimizer = "nloptwrap2",
                                               optCtrl=list(maxfun=2e4)))


schedc.glmer.fil.vacancies.expertise <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                          total - perm_excep_schedc - temp_excep_schedc) ~ 
                                      pres_switch + partyswitch + start_second_term + looming_election + 
                                      sameagency + z_pres_fildist * divided + z_polarization +
                                      z_rollcalls +  prev_exec_pct + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                                      zmecompmedian + (1|agency) + (1|pres_name),
                                    data = subset(apptypes_merged, agency != "All"), family = binomial,
                                    control=glmerControl(optimizer = "nloptwrap2",
                                                         optCtrl=list(maxfun=2e4)))


schedc.glmer.fil.expertise <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                          total - perm_excep_schedc - temp_excep_schedc) ~ 
                                      pres_switch + partyswitch + start_second_term + looming_election + 
                                      sameagency + z_pres_fildist * divided + z_polarization +
                                      z_rollcalls +  prev_exec_pct + central_staff_agency + zmecompmedian + 
                                      (1|agency) + (1|pres_name),
                                    data = subset(apptypes_merged, agency != "All"), family = binomial,
                                    control=glmerControl(optimizer = "nloptwrap2",
                                                         optCtrl=list(maxfun=2e4)))

schedc.glmer.fil.vacancies.meandur <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                          total - perm_excep_schedc - temp_excep_schedc) ~ 
                                            pres_switch + partyswitch + start_second_term + looming_election + 
                                            sameagency + z_pres_fildist * divided + z_polarization +
                                      z_rollcalls + z_meanduration + prev_exec_pct + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                                      (1|agency) + (1|pres_name),
                                    data = subset(apptypes_merged, agency != "All"), family = binomial,
                                    control=glmerControl(optimizer = "nloptwrap2",
                                                         optCtrl=list(maxfun=2e4)))


schedc.glmer.fil.vacancies.expertise.meandur <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                                  total - perm_excep_schedc - temp_excep_schedc) ~ 
                                              pres_switch + partyswitch + start_second_term + looming_election + 
                                              sameagency + z_pres_fildist * divided + z_polarization +
                                              z_rollcalls + z_meanduration + prev_exec_pct + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                                              zmecompmedian + (1|agency) + (1|pres_name),
                                            data = subset(apptypes_merged, agency != "All"), family = binomial,
                                            control=glmerControl(optimizer = "nloptwrap2",
                                                                 optCtrl=list(maxfun=2e4)))


schedc.glmer.fil.expertise.meandur <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                          total - perm_excep_schedc - temp_excep_schedc) ~ 
                                            pres_switch + partyswitch + start_second_term + looming_election + 
                                            sameagency + z_pres_fildist * divided + z_polarization +
                                      z_rollcalls +  z_meanduration + prev_exec_pct + central_staff_agency + zmecompmedian + 
                                      (1|agency) + (1|pres_name),
                                    data = subset(apptypes_merged, agency != "All"), family = binomial,
                                    control=glmerControl(optimizer = "nloptwrap2",
                                                         optCtrl=list(maxfun=2e4)))


schedc.glmer.fil.def <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                total - perm_excep_schedc - temp_excep_schedc) ~ 
                                pres_switch + partyswitch + start_second_term + looming_election + 
                                sameagency + z_pres_fildist * divided + z_polarization +
                            z_rollcalls + prev_exec_pct + (1|agency) + (1|pres_name),
                          data = subset(apptypes_merged, agency %in% c("DD", "AR", "AF", "NV")), family = binomial,
                          control=glmerControl(optimizer = "nloptwrap2",
                                               optCtrl=list(maxfun=2e4)))

schedc.glm.fil.def <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                    total - perm_excep_schedc - temp_excep_schedc) ~ 
                            pres_switch + partyswitch + start_second_term + looming_election + 
                            sameagency + z_pres_fildist * divided + z_polarization +
                                z_rollcalls + prev_exec_pct + factor(agency) + factor(pres_name),
                              data = subset(apptypes_merged, agency %in% c("DD", "AR", "AF", "NV")), family = quasibinomial)


schedc.glmer.fil.civ <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                    total - perm_excep_schedc - temp_excep_schedc) ~ 
                                pres_switch + partyswitch + start_second_term + looming_election + 
                                sameagency + z_pres_fildist * divided + z_polarization +
                                z_rollcalls + prev_exec_pct + central_staff_agency + (1|agency) + (1|pres_name),
                              data = subset(apptypes_merged, !(agency %in% c("All", "DD", "AR", "AF", "NV"))), family = binomial,
                              control=glmerControl(optimizer = "nloptwrap2",
                                                   optCtrl=list(maxfun=2e4)))

schedc.glmer.fil.civ.expertise <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                    total - perm_excep_schedc - temp_excep_schedc) ~ 
                                pres_switch + partyswitch + start_second_term + looming_election + 
                                sameagency + z_pres_fildist * divided + z_polarization +
                                z_rollcalls + prev_exec_pct + central_staff_agency + zmecompmedian + (1|agency) + (1|pres_name),
                              data = subset(apptypes_merged, !(agency %in% c("All", "DD", "AR", "AF", "NV"))), family = binomial,
                              control=glmerControl(optimizer = "nloptwrap2",
                                                   optCtrl=list(maxfun=2e4)))


schedc.glmer.fil.meandur <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                 total - perm_excep_schedc - temp_excep_schedc) ~ 
                                   pres_switch + partyswitch + start_second_term + looming_election + 
                                   sameagency + z_pres_fildist * divided + z_polarization +
                             z_rollcalls + z_meanduration + prev_exec_pct + (1|agency) + (1|pres_name),
                            data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                            family = binomial,
                           control=glmerControl(optimizer = "nloptwrap2",
                                                optCtrl=list(maxfun=2e4)))



schedc.glm.fil.cj <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                               total - perm_excep_schedc - temp_excep_schedc) ~ 
                           partyswitch + start_second_term + looming_election + 
                           z_pres_agencydist + z_pres_fildist * divided + z_polarization +
                           z_rollcalls + prev_exec_pct + factor(agency) + factor(pres_name),
                         data = subset(apptypes_merged, agency != "All"), family = quasibinomial)

schedc.glm.fil.cj.expertise <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                               total - perm_excep_schedc - temp_excep_schedc) ~ 
                           partyswitch + start_second_term + looming_election + 
                           z_pres_agencydist + z_pres_fildist * divided + z_polarization +
                           z_rollcalls + prev_exec_pct + zmecompmedian + factor(agency) + factor(pres_name),
                         data = subset(apptypes_merged, agency != "All"), family = quasibinomial)


schedc.glm.fil.cj.meandur <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                            total - perm_excep_schedc - temp_excep_schedc) ~ 
                           partyswitch + start_second_term + looming_election + 
                          z_pres_agencydist + z_pres_fildist * divided + z_polarization +
                        z_rollcalls + z_meanduration + prev_exec_pct + factor(agency) + factor(pres_name),
                        data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                        family = quasibinomial)

schedc.glm.fil.cj.expertise.meandur <- glm(cbind(perm_excep_schedc + temp_excep_schedc, 
                                       total - perm_excep_schedc - temp_excep_schedc) ~ 
                                   partyswitch + start_second_term + looming_election + 
                                   z_pres_agencydist + z_pres_fildist * divided + z_polarization +
                                   z_rollcalls + z_meanduration + prev_exec_pct + zmecompmedian + factor(agency) + factor(pres_name),
                                 data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                                 family = quasibinomial)


schedc.glmer.fil.cj <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                   total - perm_excep_schedc - temp_excep_schedc) ~ 
                               partyswitch + start_second_term + looming_election + 
                               z_pres_agencydist + z_pres_fildist * divided + z_polarization +
                               z_rollcalls + prev_exec_pct + central_staff_agency + (1|agency) + (1|pres_name),
                             data = subset(apptypes_merged, agency != "All"), family = binomial,
                             control=glmerControl(optimizer = "nloptwrap2",
                                                  optCtrl=list(maxfun=2e4)))

schedc.glmer.fil.cj.expertise <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                   total - perm_excep_schedc - temp_excep_schedc) ~ 
                               partyswitch + start_second_term + looming_election + 
                               z_pres_agencydist + z_pres_fildist * divided + z_polarization +
                               z_rollcalls + prev_exec_pct + central_staff_agency + zmecompmedian + (1|agency) + (1|pres_name),
                             data = subset(apptypes_merged, agency != "All"), family = binomial,
                             control=glmerControl(optimizer = "nloptwrap2",
                                                  optCtrl=list(maxfun=2e4)))


schedc.glmer.fil.cj.meandur <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                    total - perm_excep_schedc - temp_excep_schedc) ~ 
                               partyswitch + start_second_term + looming_election + 
                                z_pres_agencydist + z_pres_fildist * divided + z_polarization +
                                z_rollcalls + z_meanduration + prev_exec_pct + (1|agency) + (1|pres_name),
                               data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                               family = binomial,
                              control=glmerControl(optimizer = "nloptwrap2",
                                                   optCtrl=list(maxfun=2e4)))


schedc.glmer.fil.cj.expertise.meandur <- glmer(cbind(perm_excep_schedc + temp_excep_schedc, 
                                           total - perm_excep_schedc - temp_excep_schedc) ~ 
                                       partyswitch + start_second_term + looming_election + 
                                       z_pres_agencydist + z_pres_fildist * divided + z_polarization +
                                       z_rollcalls + z_meanduration + prev_exec_pct + zmecompmedian + (1|agency) + (1|pres_name),
                                     data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB"))), 
                                     family = binomial,
                                     control=glmerControl(optimizer = "nloptwrap2",
                                                          optCtrl=list(maxfun=2e4)))


ses.glm.cong.simple <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                      z_pres_congdist * divided + factor(agency) + factor(pres_name),
                    data = subset(apptypes_merged, agency != "All" & 
                                    year(date) > 1979), family = quasibinomial)


ses.glm.cong.simple.meandur <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                              z_pres_congdist * divided + z_meanduration + 
                              factor(agency) + factor(pres_name),
                            data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                            year(date) > 1979), family = quasibinomial)


ses.glm.cong.simple.nocsa <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                             z_pres_congdist * divided + factor(agency) + factor(pres_name),
                           data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS")) & 
                                           year(date) > 1979), family = quasibinomial)



ses.glm.cong <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                      pres_switch + partyswitch + start_second_term + looming_election + 
                      sameagency + z_pres_congdist * divided + z_polarization +
                      z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + factor(agency) + factor(pres_name),
                    data = subset(apptypes_merged, agency != "All" & 
                                    year(date) > 1979), family = quasibinomial)

ses.glm.cong.nocsa <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                            pres_switch + partyswitch + start_second_term + looming_election + 
                            sameagency + z_pres_congdist * divided + z_polarization +
                      z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + factor(agency) + factor(pres_name),
                    data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS")) & 
                                    year(date) > 1979), family = quasibinomial)


ses.glm.cong.vacancies <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                pres_switch + partyswitch + start_second_term + looming_election + 
                                sameagency + z_pres_congdist * divided + z_polarization +
                      z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                        factor(agency) + factor(pres_name),
                    data = subset(apptypes_merged, agency != "All" & 
                                    year(date) > 1979), family = quasibinomial)

ses.glm.cong.vacancies.expertise <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                pres_switch + partyswitch + start_second_term + looming_election + 
                                sameagency + z_pres_congdist * divided + z_polarization +
                                z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                                zmecompmedian + factor(agency) + factor(pres_name),
                              data = subset(apptypes_merged, agency != "All" & 
                                              year(date) > 1979), family = quasibinomial)


ses.glm.cong.expertise <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                pres_switch + partyswitch + start_second_term + looming_election + 
                                sameagency + z_pres_congdist * divided + z_polarization +
                                z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + zmecompmedian + 
                                factor(agency) + factor(pres_name),
                              data = subset(apptypes_merged, agency != "All" & 
                                              year(date) > 1979), family = quasibinomial)

ses.glm.cong.vacancies.meandur <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                        pres_switch + partyswitch + start_second_term + looming_election + 
                                        sameagency + z_pres_congdist * divided + z_polarization +
                                z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                                  I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                                factor(agency) + factor(pres_name),
                              data = subset(apptypes_merged, agency != "All" & 
                                              year(date) > 1979), family = quasibinomial)

ses.glm.cong.vacancies.expertise.meandur <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                        pres_switch + partyswitch + start_second_term + looming_election + 
                                        sameagency + z_pres_congdist * divided + z_polarization +
                                        z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                                        I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                                        zmecompmedian + factor(agency) + factor(pres_name),
                                      data = subset(apptypes_merged, agency != "All" & 
                                                      year(date) > 1979), family = quasibinomial)


ses.glm.cong.expertise.meandur <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                        pres_switch + partyswitch + start_second_term + looming_election + 
                                        sameagency + z_pres_congdist * divided + z_polarization +
                                z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + zmecompmedian + 
                                factor(agency) + factor(pres_name),
                              data = subset(apptypes_merged, agency != "All" & 
                                              year(date) > 1979), family = quasibinomial)


ses.glm.cong.meandur <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                              pres_switch + partyswitch + start_second_term + looming_election + 
                              sameagency + z_pres_congdist * divided + z_polarization +
                          z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                            factor(agency) + factor(pres_name),
                        data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                        year(date) > 1979), family = quasibinomial)



ses.glmer.cong.simple <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                          z_pres_congdist * divided + (1|agency) + (1|pres_name),
                        data = subset(apptypes_merged, agency != "All" & 
                                        year(date) > 1979), family = binomial,
                        control=glmerControl(optimizer = "nloptwrap2"))

ses.glmer.cong.simple.nocsa <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                 z_pres_congdist * divided + (1|agency) + (1|pres_name),
                               data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS")) & 
                                               year(date) > 1979), family = binomial,
                               control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.cong.simple.def <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                 z_pres_congdist * divided + (1|agency) + (1|pres_name),
                               data = subset(apptypes_merged, agency %in% c("DD", "AR", "AF", "NV") & 
                                               year(date) > 1979), family = binomial,
                               control=glmerControl(optimizer = "nloptwrap2"))

ses.glm.cong.simple.def <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                     z_pres_congdist * divided + factor(agency) + factor(pres_name),
                                   data = subset(apptypes_merged, agency %in% c("DD", "AR", "AF", "NV") & 
                                                   year(date) > 1979), family = quasibinomial)


ses.glmer.cong.simple.civ <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                     z_pres_congdist * divided + (1|agency) + (1|pres_name),
                                   data = subset(apptypes_merged, !(agency %in% c("All", "DD", "AR", "AF", "NV")) & 
                                                   year(date) > 1979), family = binomial,
                                   control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.cong.only.meandur <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                         z_meanduration + (1|agency) + (1|pres_name),
                                       data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                                       year(date) > 1979), family = binomial,
                                       control=glmerControl(optimizer = "nloptwrap2"))

ses.glm.cong.only.meandur <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                       z_meanduration + factor(agency) + factor(pres_name),
                                     data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                                     year(date) > 1979), family = quasibinomial)


ses.glmer.cong.simple.meandur <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                 z_pres_congdist * divided + z_meanduration + (1|agency) + (1|pres_name),
                               data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                               year(date) > 1979), family = binomial,
                               control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.cong.meandur <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                  pres_switch + partyswitch + start_second_term + looming_election + 
                                  sameagency + z_pres_congdist * divided + z_polarization +
                          z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + (1|agency) + (1|pres_name),
                        data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                        year(date) > 1979), family = binomial,
                        control=glmerControl(optimizer = "nloptwrap2"))

ses.glmer.cong <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                          pres_switch + partyswitch + start_second_term + looming_election + 
                          sameagency + z_pres_congdist * divided + z_polarization +
                          z_rollcalls + central_staff_agency + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                          (1|agency) + (1|pres_name),
                       data = subset(apptypes_merged, agency != "All" & 
                                       year(date) > 1979), family = binomial,
                       control=glmerControl(optimizer = "nloptwrap2"))

ses.glmer.cong.nocsa <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                pres_switch + partyswitch + start_second_term + looming_election + 
                                sameagency + z_pres_congdist * divided + z_polarization +
                          z_rollcalls + central_staff_agency + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                          (1|agency) + (1|pres_name),
                        data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS")) & 
                                        year(date) > 1979), family = binomial,
                        control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.cong.nocsa.expertise <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                pres_switch + partyswitch + start_second_term + looming_election + 
                                sameagency + z_pres_congdist * divided + z_polarization +
                                z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + zmecompmedian + 
                                (1|agency) + (1|pres_name),
                              data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS")) & 
                                              year(date) > 1979), family = binomial,
                              control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.cong.vacancies <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                    pres_switch + partyswitch + start_second_term + looming_election + 
                                    sameagency + z_pres_congdist * divided + z_polarization +
                          z_rollcalls + central_staff_agency + prev_exec_pct + prev_ses_noncareer_pct_gap + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                          (1|agency) + (1|pres_name),
                        data = subset(apptypes_merged, agency != "All" & 
                                        year(date) > 1979), family = binomial,
                        control=glmerControl(optimizer = "nloptwrap2"))

ses.glmer.cong.vacancies.expertise <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                    pres_switch + partyswitch + start_second_term + looming_election + 
                                    sameagency + z_pres_congdist * divided + z_polarization +
                                    z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                                    zmecompmedian + (1|agency) + (1|pres_name),
                                  data = subset(apptypes_merged, agency != "All" & 
                                                  year(date) > 1979), family = binomial,
                                  control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.cong.expertise <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                    pres_switch + partyswitch + start_second_term + looming_election + 
                                    sameagency + z_pres_congdist * divided + z_polarization +
                          z_rollcalls + central_staff_agency + prev_exec_pct + prev_ses_noncareer_pct_gap + zmecompmedian + 
                          (1|agency) + (1|pres_name),
                        data = subset(apptypes_merged, agency != "All" & 
                                        year(date) > 1979), family = binomial,
                        control=glmerControl(optimizer = "nloptwrap2"))

ses.glmer.cong.vacancies.meandur <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                            pres_switch + partyswitch + start_second_term + looming_election + 
                                            sameagency + z_pres_congdist * divided + z_polarization +
                                    z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                                    I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                                    (1|agency) + (1|pres_name),
                                  data = subset(apptypes_merged, agency != "All" & 
                                                  year(date) > 1979), family = binomial,
                                  control=glmerControl(optimizer = "nloptwrap2"))

ses.glmer.cong.vacancies.expertise.meandur <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                            pres_switch + partyswitch + start_second_term + looming_election + 
                                            sameagency + z_pres_congdist * divided + z_polarization +
                                            z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                                            I(prev_empty/(prev_empty + prev_interim + prev_permanent)) + 
                                            zmecompmedian + (1|agency) + (1|pres_name),
                                          data = subset(apptypes_merged, agency != "All" & 
                                                          year(date) > 1979), family = binomial,
                                          control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.cong.expertise.meandur <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                            pres_switch + partyswitch + start_second_term + looming_election + 
                                            sameagency + z_pres_congdist * divided + z_polarization +
                                    z_rollcalls + z_meanduration + prev_exec_pct + central_staff_agency + prev_ses_noncareer_pct_gap + zmecompmedian + 
                                    (1|agency) + (1|pres_name),
                                  data = subset(apptypes_merged, agency != "All" & 
                                                  year(date) > 1979), family = binomial,
                                  control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.cong.def <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                              pres_switch + partyswitch + start_second_term + looming_election + 
                              sameagency + z_pres_congdist * divided + z_polarization +
                          z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + (1|agency) + (1|pres_name),
                        data = subset(apptypes_merged, agency %in% c("DD", "AR", "AF", "NV") & 
                                        year(date) > 1979), family = binomial,
                        control=glmerControl(optimizer = "nloptwrap2"))

ses.glm.cong.def <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                          pres_switch + partyswitch + start_second_term + looming_election + 
                          sameagency + z_pres_congdist * divided + z_polarization +
                              z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + factor(agency) + factor(pres_name),
                            data = subset(apptypes_merged, agency %in% c("DD", "AR", "AF", "NV") & 
                                            year(date) > 1979), family = quasibinomial)


ses.glmer.cong.civ <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                              pres_switch + partyswitch + start_second_term + looming_election + 
                              sameagency + z_pres_congdist * divided + z_polarization +
                              z_rollcalls + central_staff_agency + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                              (1|agency) + (1|pres_name),
                            data = subset(apptypes_merged, !(agency %in% c("All", "DD", "AR", "AF", "NV")) & 
                                            year(date) > 1979), family = binomial,
                            control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.cong.civ.expertise <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                              pres_switch + partyswitch + start_second_term + looming_election + 
                              sameagency + z_pres_congdist * divided + z_polarization +
                              z_rollcalls + central_staff_agency + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                              zmecompmedian + (1|agency) + (1|pres_name),
                            data = subset(apptypes_merged, !(agency %in% c("All", "DD", "AR", "AF", "NV")) & 
                                            year(date) > 1979), family = binomial,
                            control=glmerControl(optimizer = "nloptwrap2"))


ses.glm.cong.cj.meandur <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                         partyswitch + start_second_term + looming_election + 
                         z_pres_agencydist + z_pres_congdist * divided + z_polarization +
                         z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                           factor(agency) + factor(pres_name),
                       data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                       year(date) > 1979), family = quasibinomial)

ses.glm.cong.cj.expertise.meandur <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                 partyswitch + start_second_term + looming_election + 
                                 z_pres_agencydist + z_pres_congdist * divided + z_polarization +
                                 z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                                 zmecompmedian + factor(agency) + factor(pres_name),
                               data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                               year(date) > 1979), family = quasibinomial)


ses.glm.cong.cj <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                         partyswitch + start_second_term + looming_election + 
                             z_pres_agencydist + z_pres_congdist * divided + z_polarization +
                             z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                         factor(agency) + factor(pres_name),
                           data = subset(apptypes_merged, agency != "All" & 
                                           year(date) > 1979), family = quasibinomial)

ses.glm.cong.cj.expertise <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                         partyswitch + start_second_term + looming_election + 
                         z_pres_agencydist + z_pres_congdist * divided + z_polarization +
                         z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + zmecompmedian + 
                         factor(agency) + factor(pres_name),
                       data = subset(apptypes_merged, agency != "All" & 
                                       year(date) > 1979), family = quasibinomial)


ses.glmer.cong.cj <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                             partyswitch + start_second_term + looming_election + 
                             z_pres_agencydist + z_pres_congdist * divided + z_polarization +
                             z_rollcalls + central_staff_agency + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                             (1|agency) + (1|pres_name),
                           data = subset(apptypes_merged, agency != "All" & 
                                           year(date) > 1979), family = binomial,
                           control=glmerControl(optimizer = "nloptwrap2"))

ses.glmer.cong.cj.expertise <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                             partyswitch + start_second_term + looming_election + 
                             z_pres_agencydist + z_pres_congdist * divided + z_polarization +
                             z_rollcalls + central_staff_agency + prev_exec_pct + prev_ses_noncareer_pct_gap + zmecompmedian + 
                             (1|agency) + (1|pres_name),
                           data = subset(apptypes_merged, agency != "All" & 
                                           year(date) > 1979), family = binomial,
                           control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.cong.cj.meandur <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                             partyswitch + start_second_term + looming_election + 
                             z_pres_agencydist + z_pres_congdist * divided + z_polarization +
                             z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                               (1|agency) + (1|pres_name),
                        data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                        year(date) > 1979), family = binomial,
                        control=glmerControl(optimizer = "nloptwrap2"))

ses.glmer.cong.cj.expertise.meandur <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                     partyswitch + start_second_term + looming_election + 
                                     z_pres_agencydist + z_pres_congdist * divided + z_polarization +
                                     z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                                     zmecompmedian + (1|agency) + (1|pres_name),
                                   data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                                   year(date) > 1979), family = binomial,
                                   control=glmerControl(optimizer = "nloptwrap2"))


ses.glm.fil.simple <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                     z_pres_fildist * divided + 
                       factor(agency) + factor(pres_name),
                   data = subset(apptypes_merged, agency != "All" & 
                                   year(date) > 1979), family = quasibinomial)

ses.glm.fil.simple.meandur <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                    z_pres_fildist * divided + z_meanduration + 
                                     factor(agency) + factor(pres_name),
                                   data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                                   year(date) > 1979), family = quasibinomial)


ses.glm.fil.simple.nocsa <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                            z_pres_fildist * divided + 
                            factor(agency) + factor(pres_name),
                          data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS")) & 
                                          year(date) > 1979), family = quasibinomial)



ses.glm.fil <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                     pres_switch + partyswitch + start_second_term + looming_election + 
                     sameagency + z_pres_fildist * divided + z_polarization +
                     z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                     factor(agency) + factor(pres_name),
                   data = subset(apptypes_merged, agency != "All" & 
                                   year(date) > 1979), family = quasibinomial)

ses.glm.fil.nocsa <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                           pres_switch + partyswitch + start_second_term + looming_election + 
                           sameagency + z_pres_fildist * divided + z_polarization +
                     z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                     factor(agency) + factor(pres_name),
                   data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS")) & 
                                   year(date) > 1979), family = quasibinomial)


ses.glm.fil.vacancies <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                               pres_switch + partyswitch + start_second_term + looming_election + 
                               sameagency + z_pres_fildist * divided + z_polarization +
                     z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) +
                     factor(agency) + factor(pres_name),
                   data = subset(apptypes_merged, agency != "All" & 
                                   year(date) > 1979), family = quasibinomial)


ses.glm.fil.vacancies.expertise <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                               pres_switch + partyswitch + start_second_term + looming_election + 
                               sameagency + z_pres_fildist * divided + z_polarization +
                               z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) +
                               zmecompmedian + factor(agency) + factor(pres_name),
                             data = subset(apptypes_merged, agency != "All" & 
                                             year(date) > 1979), family = quasibinomial)


ses.glm.fil.expertise <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                               pres_switch + partyswitch + start_second_term + looming_election + 
                               sameagency + z_pres_fildist * divided + z_polarization +
                     z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + zmecompmedian + 
                     factor(agency) + factor(pres_name),
                   data = subset(apptypes_merged, agency != "All" & 
                                   year(date) > 1979), family = quasibinomial)

ses.glm.fil.vacancies.meandur <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                       pres_switch + partyswitch + start_second_term + looming_election + 
                                       sameagency + z_pres_fildist * divided + z_polarization +
                               z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                                 I(prev_empty/(prev_empty + prev_interim + prev_permanent)) +
                               factor(agency) + factor(pres_name),
                             data = subset(apptypes_merged, agency != "All" & 
                                             year(date) > 1979), family = quasibinomial)

ses.glm.fil.vacancies.expertise.meandur <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                       pres_switch + partyswitch + start_second_term + looming_election + 
                                       sameagency + z_pres_fildist * divided + z_polarization +
                                       z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                                       I(prev_empty/(prev_empty + prev_interim + prev_permanent)) +
                                       zmecompmedian + factor(agency) + factor(pres_name),
                                     data = subset(apptypes_merged, agency != "All" & 
                                                     year(date) > 1979), family = quasibinomial)


ses.glm.fil.expertise.meandur <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                       pres_switch + partyswitch + start_second_term + looming_election + 
                                       sameagency + z_pres_fildist * divided + z_polarization +
                               z_rollcalls + z_meanduration + prev_exec_pct +prev_ses_noncareer_pct_gap + zmecompmedian + 
                               factor(agency) + factor(pres_name),
                             data = subset(apptypes_merged, agency != "All" & 
                                             year(date) > 1979), family = quasibinomial)


ses.glm.fil.meandur <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                             pres_switch + partyswitch + start_second_term + looming_election + 
                             sameagency + z_pres_fildist * divided + z_polarization +
                         z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                           factor(agency) + factor(pres_name),
                       data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                       year(date) > 1979), family = quasibinomial)


ses.glmer.fil.simple <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                        z_pres_fildist * divided + (1|agency) + (1|pres_name),
                       data = subset(apptypes_merged, agency != "All" & 
                                       year(date) > 1979), family = binomial,
                       control=glmerControl(optimizer = "nloptwrap2"))

ses.glmer.fil.simple.nocsa <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                z_pres_fildist * divided + (1|agency) + (1|pres_name),
                              data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS")) & 
                                              year(date) > 1979), family = binomial,
                              control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.fil.simple.def <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                z_pres_fildist * divided + (1|agency) + (1|pres_name),
                              data = subset(apptypes_merged, agency %in% c("DD", "AR", "AF", "NV") & 
                                              year(date) > 1979), family = binomial,
                              control=glmerControl(optimizer = "nloptwrap2"))

ses.glm.fil.simple.def <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                    z_pres_fildist * divided + factor(agency) + factor(pres_name),
                                  data = subset(apptypes_merged, agency %in% c("DD", "AR", "AF", "NV") & 
                                                  year(date) > 1979), family = quasibinomial)


ses.glmer.fil.simple.civ <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                    z_pres_fildist * divided + (1|agency) + (1|pres_name),
                                  data = subset(apptypes_merged, !(agency %in% c("All", "DD", "AR", "AF", "NV")) & 
                                                  year(date) > 1979), family = binomial,
                                  control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.fil.only.meandur <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                        z_meanduration + (1|agency) + (1|pres_name),
                                      data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                                      year(date) > 1979), family = binomial,
                                      control=glmerControl(optimizer = "nloptwrap2"))

ses.glm.fil.only.meandur <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                      z_meanduration + factor(agency) + factor(pres_name),
                                    data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                                    year(date) > 1979), family = quasibinomial)


ses.glmer.fil.simple.meandur <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                z_pres_fildist * divided + z_meanduration + (1|agency) + (1|pres_name),
                              data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                              year(date) > 1979), family = binomial,
                              control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.fil <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                         pres_switch + partyswitch + start_second_term + looming_election + 
                         sameagency + z_pres_fildist * divided + z_polarization +
                         z_rollcalls + central_staff_agency + prev_exec_pct + prev_ses_noncareer_pct_gap +
                         (1|agency) + (1|pres_name),
                       data = subset(apptypes_merged, agency != "All" & 
                                       year(date) > 1979), family = binomial,
                       control=glmerControl(optimizer = "nloptwrap2"))

ses.glmer.fil.nocsa <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                               pres_switch + partyswitch + start_second_term + looming_election + 
                               sameagency + z_pres_fildist * divided + z_polarization +
                         z_rollcalls + central_staff_agency + prev_exec_pct + prev_ses_noncareer_pct_gap +
                         (1|agency) + (1|pres_name),
                       data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS")) & 
                                       year(date) > 1979), family = binomial,
                       control=glmerControl(optimizer = "nloptwrap2"))



ses.glmer.fil.nocsa.expertise <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                               pres_switch + partyswitch + start_second_term + looming_election + 
                               sameagency + z_pres_fildist * divided + z_polarization +
                               z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + zmecompmedian + 
                               (1|agency) + (1|pres_name),
                             data = subset(apptypes_merged, !(agency %in% c("All", "BO", "OM", "GS")) & 
                                             year(date) > 1979), family = binomial,
                             control=glmerControl(optimizer = "nloptwrap2"))

ses.glmer.fil.vacancies <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                   pres_switch + partyswitch + start_second_term + looming_election + 
                                   sameagency + z_pres_fildist * divided + z_polarization +
                         z_rollcalls  + prev_exec_pct + prev_ses_noncareer_pct_gap  + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) +
                         (1|agency) + (1|pres_name),
                       data = subset(apptypes_merged, agency != "All" & 
                                       year(date) > 1979), family = binomial,
                       control=glmerControl(optimizer = "nloptwrap2"))

ses.glmer.fil.vacancies.expertise <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                   pres_switch + partyswitch + start_second_term + looming_election + 
                                   sameagency + z_pres_fildist * divided + z_polarization +
                                   z_rollcalls  + prev_exec_pct + prev_ses_noncareer_pct_gap  + I(prev_empty/(prev_empty + prev_interim + prev_permanent)) +
                                   zmecompmedian + (1|agency) + (1|pres_name),
                                 data = subset(apptypes_merged, agency != "All" & 
                                                 year(date) > 1979), family = binomial,
                                 control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.fil.expertise <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                   pres_switch + partyswitch + start_second_term + looming_election + 
                                   sameagency + z_pres_fildist * divided + z_polarization +
                         z_rollcalls + central_staff_agency + prev_exec_pct + prev_ses_noncareer_pct_gap + zmecompmedian +
                         (1|agency) + (1|pres_name),
                       data = subset(apptypes_merged, agency != "All" & 
                                       year(date) > 1979), family = binomial,
                       control=glmerControl(optimizer = "nloptwrap2"))

ses.glmer.fil.vacancies.meandur <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                           pres_switch + partyswitch + start_second_term + looming_election + 
                                           sameagency + z_pres_fildist * divided + z_polarization +
                                   z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap  + 
                                     I(prev_empty/(prev_empty + prev_interim + prev_permanent)) +
                                   (1|agency) + (1|pres_name),
                                 data = subset(apptypes_merged, agency != "All" & 
                                                 year(date) > 1979), family = binomial,
                                 control=glmerControl(optimizer = "nloptwrap2"))

ses.glmer.fil.vacancies.expertise.meandur <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                           pres_switch + partyswitch + start_second_term + looming_election + 
                                           sameagency + z_pres_fildist * divided + z_polarization +
                                           z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap  + 
                                           I(prev_empty/(prev_empty + prev_interim + prev_permanent)) +
                                           zmecompmedian + (1|agency) + (1|pres_name),
                                         data = subset(apptypes_merged, agency != "All" & 
                                                         year(date) > 1979), family = binomial,
                                         control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.fil.expertise.meandur <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                           pres_switch + partyswitch + start_second_term + looming_election + 
                                           sameagency + z_pres_fildist * divided + z_polarization +
                                   z_rollcalls + z_meanduration + prev_exec_pct + central_staff_agency + prev_ses_noncareer_pct_gap + zmecompmedian +
                                   (1|agency) + (1|pres_name),
                                 data = subset(apptypes_merged, agency != "All" & 
                                                 year(date) > 1979), family = binomial,
                                 control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.fil.def <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                             pres_switch + partyswitch + start_second_term + looming_election + 
                             sameagency + z_pres_fildist * divided + z_polarization +
                         z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + (1|agency) + (1|pres_name),
                       data = subset(apptypes_merged, agency %in% c("DD", "AR", "AF", "NV") & 
                                       year(date) > 1979), family = binomial,
                       control=glmerControl(optimizer = "nloptwrap2"))

ses.glm.fil.def <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                         pres_switch + partyswitch + start_second_term + looming_election + 
                         sameagency + z_pres_fildist * divided + z_polarization +
                             z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + factor(agency) + factor(pres_name),
                           data = subset(apptypes_merged, agency %in% c("DD", "AR", "AF", "NV") & 
                                           year(date) > 1979), family = quasibinomial)



ses.glmer.fil.civ <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                             pres_switch + partyswitch + start_second_term + looming_election + 
                             sameagency + z_pres_fildist * divided + z_polarization +
                             z_rollcalls + central_staff_agency + prev_exec_pct + prev_ses_noncareer_pct_gap +
                             (1|agency) + (1|pres_name),
                           data = subset(apptypes_merged, !(agency %in% c("All", "DD", "AR", "AF", "NV")) & 
                                           year(date) > 1979), family = binomial,
                           control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.fil.civ.expertise <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                             pres_switch + partyswitch + start_second_term + looming_election + 
                             sameagency + z_pres_fildist * divided + z_polarization +
                             z_rollcalls + central_staff_agency + prev_exec_pct + prev_ses_noncareer_pct_gap +
                             zmecompmedian + (1|agency) + (1|pres_name),
                           data = subset(apptypes_merged, !(agency %in% c("All", "DD", "AR", "AF", "NV")) & 
                                           year(date) > 1979), family = binomial,
                           control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.fil.meandur <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                 pres_switch + partyswitch + start_second_term + looming_election + 
                                 sameagency + z_pres_fildist * divided + z_polarization +
                          z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                            (1|agency) + (1|pres_name),
                        data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                        year(date) > 1979), family = binomial,
                        control=glmerControl(optimizer = "nloptwrap2"))

ses.glm.fil.cj <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                        partyswitch + start_second_term + looming_election + 
                        z_pres_agencydist + z_pres_fildist * divided + z_polarization +
                        z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                        factor(agency) + factor(pres_name),
                      data = subset(apptypes_merged, agency != "All" & 
                                      year(date) > 1979), family = quasibinomial)

ses.glm.fil.cj.expertise <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                        partyswitch + start_second_term + looming_election + 
                        z_pres_agencydist + z_pres_fildist * divided + z_polarization +
                        z_rollcalls + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                        zmecompmedian + factor(agency) + factor(pres_name),
                      data = subset(apptypes_merged, agency != "All" & 
                                      year(date) > 1979), family = quasibinomial)


ses.glm.fil.cj.meandur <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                        partyswitch + start_second_term + looming_election + 
                            z_pres_agencydist + z_pres_fildist * divided + z_polarization +
                            z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                          factor(agency) + factor(pres_name),
                          data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                          year(date) > 1979), family = quasibinomial)


ses.glm.fil.cj.expertise.meandur <- glm(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                partyswitch + start_second_term + looming_election + 
                                z_pres_agencydist + z_pres_fildist * divided + z_polarization +
                                z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                                zmecompmedian + factor(agency) + factor(pres_name),
                              data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                              year(date) > 1979), family = quasibinomial)


ses.glmer.fil.cj <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                            partyswitch + start_second_term + looming_election + 
                             z_pres_agencydist + z_pres_fildist * divided + z_polarization +
                             z_rollcalls + central_staff_agency + prev_exec_pct + prev_ses_noncareer_pct_gap +
                            (1|agency) + (1|pres_name),
                           data = subset(apptypes_merged, agency != "All" & 
                                           year(date) > 1979), family = binomial,
                           control=glmerControl(optimizer = "nloptwrap2"))

ses.glmer.fil.cj.expertise <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                            partyswitch + start_second_term + looming_election + 
                            z_pres_agencydist + z_pres_fildist * divided + z_polarization +
                            z_rollcalls + central_staff_agency + prev_exec_pct + prev_ses_noncareer_pct_gap + zmecompmedian + 
                            (1|agency) + (1|pres_name),
                          data = subset(apptypes_merged, agency != "All" & 
                                          year(date) > 1979), family = binomial,
                          control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.fil.cj.meandur <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                            partyswitch + start_second_term + looming_election + 
                            z_pres_agencydist + z_pres_fildist * divided + z_polarization +
                            z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                              (1|agency) + (1|pres_name),
                          data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                          year(date) > 1979), family = binomial,
                          control=glmerControl(optimizer = "nloptwrap2"))


ses.glmer.fil.cj.expertise.meandur <- glmer(cbind(perm_ses_noncareer, perm_ses_career) ~ 
                                    partyswitch + start_second_term + looming_election + 
                                    z_pres_agencydist + z_pres_fildist * divided + z_polarization +
                                    z_rollcalls + z_meanduration + prev_exec_pct + prev_ses_noncareer_pct_gap + 
                                    zmecompmedian + (1|agency) + (1|pres_name),
                                  data = subset(apptypes_merged, !(agency %in% c("All", "OPM", "GSA", "OMB")) & 
                                                  year(date) > 1979), family = binomial,
                                  control=glmerControl(optimizer = "nloptwrap2"))

# utility functions to extract correct test statistics for htmlreg files
source("HR-Utility-Functions.R")


# export to tables
htmlreg(file = "Table-1.doc",inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        list(schedc.glmer.cong.simple,
             schedc.glmer.cong.expertise, 
             schedc.glmer.fil.simple,
             schedc.glmer.fil.expertise, 
             ses.glmer.cong.simple, 
             ses.glmer.cong.expertise, 
             ses.glmer.fil.simple,
             ses.glmer.fil.expertise),
        custom.coef.names = c("Constant", 
                              "President-Senate Distance",
                              "Divided Government",
                              "President-Senate Distance times Divided Government",
                              "New President",
                              "New Presidential Party",
                              "Presidential Reelection",
                              "Upcoming Presidential Election",
                              "President-Agency Alignment",
                              "Congressional Polarization",
                              "Quarterly Senate Roll Calls",
                              "Previous Quarter Executive Schedule Percentage",
                              "Central Staffing Agency",
                              "Managerial Expertise",
                              "President-Filibuster Distance",
                              "President-Filibuster Distance times Divided Government",
                              "Previous Quarter SES Noncareer Percentage Gap"),
        digits = 3, stars = c(0.01, 0.05, 0.1),        
        omit.coef = c("factor"))


htmlreg(file = "Table-A1.doc",inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        list(schedc.glm.cong.simple,
             schedc.glm.cong.expertise, 
             schedc.glm.fil.simple,
             schedc.glm.fil.expertise, 
             ses.glm.cong.simple, 
             ses.glm.cong.expertise, 
             ses.glm.fil.simple,
             ses.glm.fil.expertise),
        custom.coef.names = c("Constant", 
                              "President-Senate Distance",
                              "Divided Government",
                              "President-Senate Distance times Divided Government",
                              "New President",
                              "New Presidential Party",
                              "Presidential Reelection",
                              "Upcoming Presidential Election",
                              "President-Agency Alignment",
                              "Congressional Polarization",
                              "Quarterly Senate Roll Calls",
                              "Previous Quarter Executive Schedule Percentage",
                              "Managerial Expertise",
                              "President-Filibuster Distance",
                              "President-Filibuster Distance times Divided Government",
                              "Previous Quarter SES Noncareer Percentage Gap"),
        digits = 3, stars = c(0.01, 0.05, 0.1),        
        omit.coef = c("factor"))



htmlreg(file = "Table-A2.doc",inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        list(schedc.glmer.cong.simple.meandur,
             schedc.glmer.cong.expertise.meandur, 
             schedc.glmer.fil.simple.meandur,
             schedc.glmer.fil.expertise.meandur, 
             ses.glmer.cong.simple.meandur, 
             ses.glmer.cong.expertise.meandur, 
             ses.glmer.fil.simple.meandur,
             ses.glmer.fil.expertise.meandur),
        custom.coef.names = c("Constant", 
                              "President-Senate Distance",
                              "Divided Government",
                              "Average Confirmation Duration",
                              "President-Senate Distance times Divided Government",
                              "New President",
                              "New Presidential Party",
                              "Presidential Reelection",
                              "Upcoming Presidential Election",
                              "President-Agency Alignment",
                              "Congressional Polarization",
                              "Quarterly Senate Roll Calls",
                              "Previous Quarter Executive Schedule Percentage",
                              "Managerial Expertise",
                              "President-Filibuster Distance",
                              "President-Filibuster Distance times Divided Government",
                              "Previous Quarter SES Noncareer Percentage Gap"),
        digits = 3, stars = c(0.01, 0.05, 0.1),        
        omit.coef = c("factor"))


htmlreg(file = "Table-A3.doc",inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        list(schedc.glm.cong.simple.meandur,
             schedc.glm.cong.expertise.meandur, 
             schedc.glm.fil.simple.meandur,
             schedc.glm.fil.expertise.meandur, 
             ses.glm.cong.simple.meandur, 
             ses.glm.cong.expertise.meandur, 
             ses.glm.fil.simple.meandur,
             ses.glm.fil.expertise.meandur),
        custom.coef.names = c("Constant", 
                              "President-Senate Distance",
                              "Divided Government",
                              "Average Confirmation Duration",
                              "President-Senate Distance times Divided Government",
                              "New President",
                              "New Presidential Party",
                              "Presidential Reelection",
                              "Upcoming Presidential Election",
                              "President-Agency Alignment",
                              "Congressional Polarization",
                              "Quarterly Senate Roll Calls",
                              "Previous Quarter Executive Schedule Percentage",
                              "Managerial Expertise",
                              "President-Filibuster Distance",
                              "President-Filibuster Distance times Divided Government",
                              "Previous Quarter SES Noncareer Percentage Gap"),
        digits = 3, stars = c(0.01, 0.05, 0.1),        
        omit.coef = c("factor"))


htmlreg(file = "Table-A4.doc",inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        list(schedc.glmer.cong.cj.expertise, 
             schedc.glmer.cong.cj.expertise.meandur,
             schedc.glmer.fil.cj.expertise, 
             schedc.glmer.fil.cj.expertise.meandur, 
             ses.glmer.cong.cj.expertise, 
             ses.glmer.cong.cj.expertise.meandur, 
             ses.glmer.fil.cj.expertise,
             ses.glmer.fil.cj.expertise.meandur),
        custom.coef.names = c("Constant", 
                              "New Presidential Party",
                              "Presidential Reelection",
                              "Upcoming Presidential Election",
                              "President-Agency Distance",
                              "President-Senate Distance",
                              "Divided Government",
                              "Congressional Polarization",
                              "Quarterly Senate Roll Calls",
                              "Central Staffing Agency",
                              "Previous Quarter Executive Schedule Percentage",
                              "Managerial Expertise",
                              "President-Senate Distance times Divided Government",
                              "Average Confirmation Duration",
                              "President-Filibuster Distance",
                              "President-Filibuster Distance times Divided Government",
                              "Previous Quarter SES Noncareer Percentage Gap"),
        digits = 3, stars = c(0.01, 0.05, 0.1),        
        omit.coef = "factor")



htmlreg(file = "Table-A5.doc",inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        list(schedc.glm.cong.cj.expertise, 
             schedc.glm.cong.cj.expertise.meandur,
             schedc.glm.fil.cj.expertise, 
             schedc.glm.fil.cj.expertise.meandur, 
             ses.glm.cong.cj.expertise, 
             ses.glm.cong.cj.expertise.meandur, 
             ses.glm.fil.cj.expertise,
             ses.glm.fil.cj.expertise.meandur),
        custom.coef.names = c("Constant", 
                              "New Presidential Party",
                              "Presidential Reelection",
                              "Upcoming Presidential Election",
                              "President-Agency Distance",
                              "President-Senate Distance",
                              "Divided Government",
                              "Congressional Polarization",
                              "Quarterly Senate Roll Calls",
                              "Previous Quarter Executive Schedule Percentage",
                              "Managerial Expertise",
                              "President-Senate Distance times Divided Government",
                              "Average Confirmation Duration",
                              "President-Filibuster Distance",
                              "President-Filibuster Distance times Divided Government",
                              "Previous Quarter SES Noncareer Percentage Gap"),
        digits = 3, stars = c(0.01, 0.05, 0.1),        
        omit.coef = "factor")


htmlreg(file = "Table-A6.doc",inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        list(schedc.glmer.cong, 
             schedc.glmer.cong.meandur,
             schedc.glmer.fil, 
             schedc.glmer.fil.meandur, 
             ses.glmer.cong, 
             ses.glmer.cong.meandur,
             ses.glmer.fil, 
             ses.glmer.fil.meandur),
        custom.coef.names = c("Constant", 
                              "New President",
                              "New Presidential Party",
                              "Presidential Reelection",
                              "Upcoming Presidential Election",
                              "President-Agency Alignment",
                              "President-Senate Distance",
                              "Divided Government",
                              "Congressional Polarization",
                              "Quarterly Senate Roll Calls",
                              "Previous Quarter Executive Schedule Percentage",
                              "Central Staffing Agency",
                              "President-Senate Distance times Divided Government",
                              "Average Confirmation Duration",
                              "President-Filibuster Distance",
                              "President-Filibuster Distance times Divided Government",
                              "Previous Quarter SES Noncareer Percentage Gap"),
        digits = 3, stars = c(0.01, 0.05, 0.1),        
        omit.coefs = "factor")



htmlreg(file = "Table-A7.doc",inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        list(schedc.glm.cong, 
             schedc.glm.cong.meandur,
             schedc.glm.fil, 
             schedc.glm.fil.meandur, 
             ses.glm.cong, 
             ses.glm.cong.meandur,
             ses.glm.fil, 
             ses.glm.fil.meandur),
        custom.coef.names = c("Constant", 
                              "New President",
                              "New Presidential Party",
                              "Presidential Reelection",
                              "Upcoming Presidential Election",
                              "President-Agency Alignment",
                              "President-Senate Distance",
                              "Divided Government",
                              "Congressional Polarization",
                              "Quarterly Senate Roll Calls",
                              "Previous Quarter Executive Schedule Percentage",
                              "President-Senate Distance times Divided Government",
                              "Average Confirmation Duration",
                              "President-Filibuster Distance",
                              "President-Filibuster Distance times Divided Government",
                              "Previous Quarter SES Noncareer Percentage Gap"),
        digits = 3, stars = c(0.01, 0.05, 0.1),        
        omit.coef = "factor")


htmlreg(file = "Table-A8.doc",inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        list(schedc.glmer.cong.vacancies.expertise, 
             schedc.glmer.cong.vacancies.expertise.meandur,
             schedc.glmer.fil.vacancies.expertise, 
             schedc.glmer.fil.vacancies.expertise.meandur, 
             ses.glmer.cong.vacancies.expertise, 
             ses.glmer.cong.vacancies.expertise.meandur,
             ses.glmer.fil.vacancies.expertise, 
             ses.glmer.fil.vacancies.expertise.meandur),
        custom.coef.names = c("Constant", 
                              "New President",
                              "New Presidential Party",
                              "Presidential Reelection",
                              "Upcoming Presidential Election",
                              "President-Agency Alignment",
                              "President-Senate Distance",
                              "Divided Government",
                              "Congressional Polarization",
                              "Quarterly Senate Roll Calls",
                              "Previous Quarter Executive Schedule Percentage",
                              "Previous Year Vacancy Percentage",
                              "Managerial Expertise",
                              "President-Senate Distance times Divided Government",
                              "Average Confirmation Duration",
                              "President-Filibuster Distance",
                              "President-Filibuster Distance times Divided Government",
                              "Previous Quarter SES Noncareer Percentage Gap"),
        digits = 3, stars = c(0.01, 0.05, 0.1),        
        omit.coef = "factor")
        

htmlreg(file = "Table-A9.doc",inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        list(schedc.glm.cong.vacancies.expertise, 
             schedc.glm.cong.vacancies.expertise.meandur,
             schedc.glm.fil.vacancies.expertise, 
             schedc.glm.fil.vacancies.expertise.meandur, 
             ses.glm.cong.vacancies.expertise, 
             ses.glm.cong.vacancies.expertise.meandur,
             ses.glm.fil.vacancies.expertise, 
             ses.glm.fil.vacancies.expertise.meandur),
        custom.coef.names = c("Constant", 
                              "New President",
                              "New Presidential Party",
                              "Presidential Reelection",
                              "Upcoming Presidential Election",
                              "President-Agency Alignment",
                              "President-Senate Distance",
                              "Divided Government",
                              "Congressional Polarization",
                              "Quarterly Senate Roll Calls",
                              "Previous Quarter Executive Schedule Percentage",
                              "Previous Year Vacancy Percentage",
                              "Managerial Expertise",
                              "President-Senate Distance times Divided Government",
                              "Average Confirmation Duration",
                              "President-Filibuster Distance",
                              "President-Filibuster Distance times Divided Government",
                              "Previous Quarter SES Noncareer Percentage Gap"),
        digits = 3, stars = c(0.01, 0.05, 0.1),        
        omit.coef = "factor")


htmlreg(file = "Table-A10.doc",inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        list(schedc.glmer.cong.simple.def,
             schedc.glmer.cong.def, 
             schedc.glmer.fil.simple.def,
             schedc.glmer.fil.def, 
             ses.glmer.cong.simple.def, 
             ses.glmer.cong.def, 
             ses.glmer.fil.simple.def,
             ses.glmer.fil.def),
        custom.coef.names = c("Constant", 
                              "President-Senate Distance",
                              "Divided Government",
                              "President-Senate Distance times Divided Government",
                              "New President",
                              "New Presidential Party",
                              "Presidential Reelection",
                              "Upcoming Presidential Election",
                              "President-Agency Alignment",
                              "Congressional Polarization",
                              "Quarterly Senate Roll Calls",
                              "Previous Quarter Executive Schedule Percentage",
                              "President-Filibuster Distance",
                              "President-Filibuster Distance times Divided Government",
                              "Previous Quarter SES Noncareer Percentage Gap"),
        digits = 3, stars = c(0.01, 0.05, 0.1),        
        omit.coef = c("factor"))


htmlreg(file = "Table-A11.doc",inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        list(schedc.glmer.cong.simple.civ,
             schedc.glmer.cong.civ.expertise, 
             schedc.glmer.fil.simple.civ,
             schedc.glmer.fil.civ.expertise, 
             ses.glmer.cong.simple.civ, 
             ses.glmer.cong.civ.expertise, 
             ses.glmer.fil.simple.civ,
             ses.glmer.fil.civ.expertise),
        custom.coef.names = c("Constant", 
                              "President-Senate Distance",
                              "Divided Government",
                              "President-Senate Distance times Divided Government",
                              "New President",
                              "New Presidential Party",
                              "Presidential Reelection",
                              "Upcoming Presidential Election",
                              "President-Agency Alignment",
                              "Congressional Polarization",
                              "Quarterly Senate Roll Calls",
                              "Previous Quarter Executive Schedule Percentage",
                              "Central Staffing Agency",
                              "Managerial Expertise",
                              "President-Filibuster Distance",
                              "President-Filibuster Distance times Divided Government",
                              "Previous Quarter SES Noncareer Percentage Gap"),
        digits = 3, stars = c(0.01, 0.05, 0.1),        
        omit.coef = c("factor"))


htmlreg(file = "Table-A12.doc",inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        list(schedc.glmer.cong.simple.nocsa,
             schedc.glmer.cong.nocsa.expertise, 
             schedc.glmer.fil.simple.nocsa,
             schedc.glmer.fil.nocsa.expertise, 
             ses.glmer.cong.simple.nocsa, 
             ses.glmer.cong.nocsa.expertise, 
             ses.glmer.fil.simple.nocsa,
             ses.glmer.fil.nocsa.expertise),
        custom.coef.names = c("Constant", 
                              "President-Senate Distance",
                              "Divided Government",
                              "President-Senate Distance times Divided Government",
                              "New President",
                              "New Presidential Party",
                              "Presidential Reelection",
                              "Upcoming Presidential Election",
                              "President-Agency Alignment",
                              "Congressional Polarization",
                              "Quarterly Senate Roll Calls",
                              "Previous Quarter Executive Schedule Percentage",
                              "Managerial Expertise",
                              "President-Filibuster Distance",
                              "President-Filibuster Distance times Divided Government",
                              "Previous Quarter SES Noncareer Percentage Gap"),
        digits = 3, stars = c(0.01, 0.05, 0.1),        
        omit.coef = c("factor"))


htmlreg(file = "Table-A13.doc",inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        list(schedc.glmer.cong12.simple,
             schedc.glmer.cong12.expertise, 
             schedc.glmer.fil12.simple,
             schedc.glmer.fil12.expertise, 
             schedc.glm.cong12.simple,
             schedc.glm.cong12.expertise, 
             schedc.glm.fil12.simple,
             schedc.glm.fil12.expertise),
        custom.coef.names = c("Constant", 
                              "President-Senate Distance",
                              "Divided Government",
                              "President-Senate Distance times Divided Government",
                              "New President",
                              "New Presidential Party",
                              "Presidential Reelection",
                              "Upcoming Presidential Election",
                              "President-Agency Alignment",
                              "Congressional Polarization",
                              "Quarterly Senate Roll Calls",
                              "Previous Quarter Executive Schedule Percentage",
                              "Central Staffing Agency",
                              "Managerial Expertise",
                              "President-Filibuster Distance",
                              "President-Filibuster Distance times Divided Government"),
        digits = 3, stars = c(0.01, 0.05, 0.1),        
        omit.coef = "factor")

htmlreg(file = "Table-A14.doc",inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        list(schedc.glm.cong.simple.nodhs.nohhs, 
             schedc.glm.cong.nodhs.nohhs.expertise, 
             schedc.glm.fil.simple.nodhs.nohhs, 
             schedc.glm.fil.nodhs.nohhs.expertise, 
             schedc.glmer.cong.simple.nodhs.nohhs, 
             schedc.glmer.cong.nodhs.nohhs.expertise, 
             schedc.glmer.fil.simple.nodhs.nohhs, 
             schedc.glmer.fil.nodhs.nohhs.expertise),
        custom.coef.names = c("Constant", 
                              "President-Senate Distance",
                              "Divided Government",
                              "President-Senate Distance times Divided Government",
                              "New President",
                              "New Presidential Party",
                              "Presidential Reelection",
                              "Upcoming Presidential Election",
                              "President-Agency Alignment",
                              "Congressional Polarization",
                              "Quarterly Senate Roll Calls", 
                              "Previous Quarter Executive Schedule Percentage",
                              "Managerial Expertise",
                              "President-Filibuster Distance",
                              "President-Filibuster Distance times Divided Government", 
                              "Central Staffing Agency"),
        digits = 3, stars = c(0.01, 0.05, 0.1),        
        omit.coef = "factor")


# generate plots according to hanmer-kalkan approach
ses.glmer.fil.preds.unified <- ses.glmer.cong.preds.unified <- schedc.glmer.fil.preds.unified <- schedc.glmer.cong.preds.unified <- 
  ses.glmer.fil.preds.divided <- ses.glmer.cong.preds.divided <- schedc.glmer.fil.preds.divided <- schedc.glmer.cong.preds.divided <- 
  ses.glmer.fil.def.preds.unified <- ses.glmer.cong.def.preds.unified <- 
  schedc.glmer.fil.def.preds.unified <- schedc.glmer.cong.def.preds.unified <- 
  ses.glmer.fil.def.preds.divided <- ses.glmer.cong.def.preds.divided <- 
  schedc.glmer.fil.def.preds.divided <- schedc.glmer.cong.def.preds.divided <- 
  ses.glmer.fil.civ.preds.unified <- ses.glmer.cong.civ.preds.unified <- 
  schedc.glmer.fil.civ.preds.unified <- schedc.glmer.cong.civ.preds.unified <- 
  ses.glmer.fil.civ.preds.divided <- ses.glmer.cong.civ.preds.divided <- 
  schedc.glmer.fil.civ.preds.divided <- schedc.glmer.cong.civ.preds.divided <- 
  data.frame(Distance = seq(-2,2,length.out = 101),
                                  low95 = NA,
                                  low90 = NA,
                                  pred = NA,
                                  high90 = NA,
                                  high95 = NA)

dist.range <- seq(-2, 2, length.out = 101)



ses.fil.preds.unified <- ses.cong.preds.unified <- schedc.fil.preds.unified <- schedc.cong.preds.unified <- 
  ses.fil.preds.divided <- ses.cong.preds.divided <- schedc.fil.preds.divided <- schedc.cong.preds.divided <- 
  ses.fil.def.preds.unified <- ses.cong.def.preds.unified <- schedc.fil.def.preds.unified <- schedc.cong.def.preds.unified <- 
  ses.fil.def.preds.divided <- ses.cong.def.preds.divided <- schedc.fil.def.preds.divided <- schedc.cong.def.preds.divided <- 
  ses.fil.civ.preds.unified <- ses.cong.civ.preds.unified <- schedc.fil.civ.preds.unified <- schedc.cong.civ.preds.unified <- 
  ses.fil.civ.preds.divided <- ses.cong.civ.preds.divided <- schedc.fil.civ.preds.divided <- schedc.cong.civ.preds.divided <- NULL

ses.manage.fil.preds.unified <- ses.manage.cong.preds.unified <- schedc.manage.fil.preds.unified <- 
  schedc.manage.cong.preds.unified <- ses.manage.fil.preds.divided <- ses.manage.cong.preds.divided <- 
  schedc.manage.fil.preds.divided <- schedc.manage.cong.preds.divided <- NULL
  

for(i in 1:length(ses.glmer.fil.preds.unified$Distance)){
  foo.newdata.unified.cong <- foo.newdata.unified.fil <-  
    foo.newdata.divided.cong <- foo.newdata.divided.fil <- subset(apptypes_merged, agency != "All",
                                                                  select = c(pres_switch, partyswitch,
                                                                             start_second_term, sameagency,
                                                                             z_pres_congdist, z_polarization,
                                                                             z_pres_agencydist, z_pres_fildist,
                                                                             divided, z_rollcalls, z_appts_total, central_staff_agency,
                                                                             prev_exec_pct, prev_ses_noncareer_pct_gap, zmecompmedian,
                                                                             looming_election))

  foo.newdata.unified.fil$z_pres_fildist <- dist.range[i]
  foo.newdata.unified.cong$z_pres_congdist <- dist.range[i]
  foo.newdata.divided.fil$z_pres_fildist <- dist.range[i]
  foo.newdata.divided.cong$z_pres_congdist <- dist.range[i]
  foo.newdata.unified.fil$divided <- 0
  foo.newdata.divided.fil$divided <- 1
  foo.newdata.unified.cong$divided <- 0
  foo.newdata.divided.cong$divided <- 1
  
  foo.ses.fil.preds.unified <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.unified.fil, re.form = NA))
  foo.ses.cong.preds.unified <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.unified.cong, re.form = NA))
  foo.schedc.fil.preds.unified <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.unified.fil, re.form = NA))
  foo.schedc.cong.preds.unified <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.unified.cong, re.form = NA))
  foo.ses.fil.preds.divided <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.divided.fil, re.form = NA))
  foo.ses.cong.preds.divided <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.divided.cong, re.form = NA))
  foo.schedc.fil.preds.divided <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.divided.fil, re.form = NA))
  foo.schedc.cong.preds.divided <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.divided.cong, re.form = NA))
  
  ses.fil.preds.unified <- rbind(ses.fil.preds.unified, foo.ses.fil.preds.unified)
  ses.cong.preds.unified <- rbind(ses.cong.preds.unified, foo.ses.cong.preds.unified)
  schedc.fil.preds.unified <- rbind(schedc.fil.preds.unified, foo.schedc.fil.preds.unified)
  schedc.cong.preds.unified <- rbind(schedc.cong.preds.unified, foo.schedc.cong.preds.unified)
  ses.fil.preds.divided <- rbind(ses.fil.preds.divided, foo.ses.fil.preds.divided)
  ses.cong.preds.divided <- rbind(ses.cong.preds.divided, foo.ses.cong.preds.divided)
  schedc.fil.preds.divided <- rbind(schedc.fil.preds.divided, foo.schedc.fil.preds.divided)
  schedc.cong.preds.divided <- rbind(schedc.cong.preds.divided, foo.schedc.cong.preds.divided)

  }

# extract actual predictions
preds <- rbind(data.frame(SD = c("-2 to -1", "-1 to 0", "0 to 1", "1 to 2"),
                          t(apply(diff(ses.cong.preds.divided[c(1,26,51,76,101),]),1,
                                  quantile,c(0.025, 0.05, 0.5, 0.95, 0.975),na.rm = TRUE)),
                          Type = "Noncareer SES Executives\n(As Percentage of Total SES Executives)",
                          Ideology = "President-Senate Distance",
                          Govt = "Divided Government"),
               data.frame(SD = c("-2 to -1", "-1 to 0", "0 to 1", "1 to 2"),
                          t(apply(diff(ses.cong.preds.unified[c(1,26,51,76,101),]),1,
                                  quantile,c(0.025, 0.05, 0.5, 0.95, 0.975),na.rm = TRUE)),
                          Type = "Noncareer SES Executives\n(As Percentage of Total SES Executives)",
                          Ideology = "President-Senate Distance",
                          Govt = "Unified Government"),
               data.frame(SD = c("-2 to -1", "-1 to 0", "0 to 1", "1 to 2"),
                          t(apply(diff(ses.fil.preds.divided[c(1,26,51,76,101),]),1,
                                  quantile,c(0.025, 0.05, 0.5, 0.95, 0.975),na.rm = TRUE)),
                          Type = "Noncareer SES Executives\n(As Percentage of Total SES Executives)",
                          Ideology = "President-Filibuster Distance",
                          Govt = "Divided Government"),
               data.frame(SD = c("-2 to -1", "-1 to 0", "0 to 1", "1 to 2"),
                          t(apply(diff(ses.fil.preds.unified[c(1,26,51,76,101),]),1,
                                  quantile,c(0.025, 0.05, 0.5, 0.95, 0.975),na.rm = TRUE)),
                          Type = "Noncareer SES Executives\n(As Percentage of Total SES Executives)",
                          Ideology = "President-Filibuster Distance",
                          Govt = "Unified Government"),
               data.frame(SD = c("-2 to -1", "-1 to 0", "0 to 1", "1 to 2"),
                          t(apply(diff(schedc.cong.preds.divided[c(1,26,51,76,101),]),1,
                                  quantile,c(0.025, 0.05, 0.5, 0.95, 0.975),na.rm = TRUE)),
                          Type = "Schedule C Positions\n(As Percentage of Total Positions)",
                          Ideology = "President-Senate Distance", 
                          Govt = "Divided Government"),
               data.frame(SD = c("-2 to -1", "-1 to 0", "0 to 1", "1 to 2"),
                          t(apply(diff(schedc.cong.preds.unified[c(1,26,51,76,101),]),1,
                                  quantile,c(0.025, 0.05, 0.5, 0.95, 0.975),na.rm = TRUE)),
                          Type = "Schedule C Positions\n(As Percentage of Total Positions)",
                          Ideology = "President-Senate Distance",
                          Govt = "Unified Government"),
               data.frame(SD = c("-2 to -1", "-1 to 0", "0 to 1", "1 to 2"),
                          t(apply(diff(schedc.fil.preds.divided[c(1,26,51,76,101),]),1,
                                  quantile,c(0.025, 0.05, 0.5, 0.95, 0.975),na.rm = TRUE)),
                          Type = "Schedule C Positions\n(As Percentage of Total Positions)",
                          Ideology = "President-Filibuster Distance",
                          Govt = "Divided Government"),
               data.frame(SD = c("-2 to -1", "-1 to 0", "0 to 1", "1 to 2"),
                          t(apply(diff(schedc.fil.preds.unified[c(1,26,51,76,101),]),1,
                                  quantile,c(0.025, 0.05, 0.5, 0.95, 0.975),na.rm = TRUE)),
                          Type = "Schedule C Positions\n(As Percentage of Total Positions)",
                          Ideology = "President-Filibuster Distance",
                          Govt = "Unified Government"))

# relabel as factor for ggplot
preds$SD <- factor(preds$SD, levels = c("-2 to -1", "-1 to 0", "0 to 1", "1 to 2"))


fig2 <- ggplot(preds, aes(x = SD, y = X50., shape = Ideology, color = Ideology)) + 
                geom_point(position = position_dodge(width = 0.25)) + 
                facet_grid(Type~Govt, scales = "free_y") + 
                geom_errorbar(aes(ymin = X5., ymax = X95.), width = 0, size = 0.75, position = position_dodge(width = 0.25)) + 
                geom_errorbar(aes(ymin = X2.5., ymax = X97.5.), width = 0, position = position_dodge(width = 0.25)) + 
                geom_hline(yintercept = 0, linetype = "dotted") + 
                theme_bw() + 
                ylab('Percentage Point Change in Proportion of Indicated Position Type\n') + 
                theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
                xlab('\nStandard Deviation Change in Indicated Ideological Measure\n(Relative to Mean)') + 
                scale_shape_discrete(name = "Ideological Measure") + 
                scale_y_continuous(trans = "asinh") + 
                scale_color_grey(name = "Ideological Measure", start = 0.2, end = 0.6)

ggsave(fig2, file = "Figure-2.pdf", height = 6.25, width = 10)
ggsave(fig2, file = "Figure-2.png", height = 6.25, width = 10)



# effect of divided government
divided.preds <- rbind(data.frame(Distance = dist.range,
                                  t(apply(ses.cong.preds.divided - ses.cong.preds.unified,1,quantile,
                                          c(0.025,0.05,0.5,0.95,0.975), na.rm = TRUE)),
                                  Type = "Noncareer SES Executives\n(As Percentage of Total SES Executives)",
                                  Ideology = "President-Senate Distance",
                                  Xleft = range(apptypes_merged$z_pres_congdist[apptypes_merged$divided == 1], na.rm = TRUE)[1],
                                  Xright = range(apptypes_merged$z_pres_congdist[apptypes_merged$divided == 0], na.rm = TRUE)[2]),
                       data.frame(Distance = dist.range,
                                  t(apply(ses.fil.preds.divided - ses.fil.preds.unified,1,quantile,
                                          c(0.025,0.05,0.5,0.95,0.975), na.rm = TRUE)),
                                  Type = "Noncareer SES Executives\n(As Percentage of Total SES Executives)",
                                  Ideology = "President-Filibuster Distance",
                                  Xleft = range(apptypes_merged$z_pres_fildist[apptypes_merged$divided == 1], na.rm = TRUE)[1],
                                  Xright = range(apptypes_merged$z_pres_fildist[apptypes_merged$divided == 0], na.rm = TRUE)[2]),
                       data.frame(Distance = dist.range,
                                  t(apply(schedc.cong.preds.divided - schedc.cong.preds.unified,1,quantile,
                                          c(0.025,0.05,0.5,0.95,0.975), na.rm = TRUE)),
                                  Type = "Schedule C Positions\n(As Percentage of Total Positions)",
                                  Ideology = "President-Senate Distance",
                                  Xleft = range(apptypes_merged$z_pres_congdist[apptypes_merged$divided == 1], na.rm = TRUE)[1],
                                  Xright = range(apptypes_merged$z_pres_congdist[apptypes_merged$divided == 0], na.rm = TRUE)[2]),
                       data.frame(Distance = dist.range,
                                  t(apply(schedc.fil.preds.divided - schedc.fil.preds.unified,1,quantile,
                                          c(0.025,0.05,0.5,0.95,0.975), na.rm = TRUE)),
                                  Type = "Schedule C Positions\n(As Percentage of Total Positions)",
                                  Ideology = "President-Filibuster Distance",
                                  Xleft = range(apptypes_merged$z_pres_fildist[apptypes_merged$divided == 1], na.rm = TRUE)[1],
                                  Xright = range(apptypes_merged$z_pres_fildist[apptypes_merged$divided == 0], na.rm = TRUE)[2]))
         
# convert to factor for ggplot    
divided.preds$Ideology <- factor(divided.preds$Ideology, levels = c("President-Senate Distance",
                                                                    "President-Filibuster Distance"))

fig4 <- ggplot(divided.preds, aes(x = Distance, y = X50.)) + 
  geom_line() + 
  geom_vline(aes(xintercept = Xleft), linetype = "dashed") + 
  geom_vline(aes(xintercept = Xright), linetype = "dashed") + 
  facet_grid(Type~Ideology, scales = "free_y") + 
  geom_ribbon(aes(ymin = X2.5., ymax = X97.5.), alpha = I(0.2)) + 
  geom_ribbon(aes(ymin = X5., ymax = X95.), alpha = I(0.6)) + 
  geom_hline(yintercept = 0, linetype = "dotted") + 
  theme_bw() + 
  ylab('Percentage Point Change in Proportion of Indicated Position Type\n') + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  xlab('\nIdeology\n(Number of Standard Deviations Relative to Mean)') + 
  scale_shape_discrete(name = "Ideological Measure") + 
  scale_y_continuous(trans = "asinh") + 
  scale_color_grey(name = "Ideological Measure", start = 0.2, end = 0.6)

ggsave(fig4, file = "Figure-4.pdf", height = 6, width = 10)
ggsave(fig4, file = "Figure-4.png", height = 6, width = 10)


# predictions in terms of number of appointees (figure 3)
ses.glmer.fil.preds.unified.total <- 
  ses.glmer.cong.preds.unified.total <- 
  schedc.glmer.fil.preds.unified.total <- 
  schedc.glmer.cong.preds.unified.total <- 
  ses.glmer.fil.preds.divided.total <- 
  ses.glmer.cong.preds.divided.total <- 
  schedc.glmer.fil.preds.divided.total <- 
  schedc.glmer.cong.preds.divided.total <- 
  data.frame(Distance = seq(-2,2,length.out = 101),
             low95 = NA,
             low90 = NA,
             pred = NA,
             high90 = NA,
             high95 = NA)

dist.range <- seq(-2, 2, length.out = 101)

ses.fil.preds.unified.total <- ses.cong.preds.unified.total <- 
  schedc.fil.preds.unified.total <- schedc.cong.preds.unified.total <- 
  ses.fil.preds.divided.total <- ses.cong.preds.divided.total <- 
  schedc.fil.preds.divided.total <- schedc.cong.preds.divided.total <- NULL

for(i in 1:length(ses.glmer.fil.preds.unified.total$Distance)){
  foo.newdata.unified.cong.total <- foo.newdata.unified.fil.total <-  
    foo.newdata.divided.cong.total <- foo.newdata.divided.fil.total <- subset(apptypes_merged, agency != "All",
                                                                  select = c(pres_switch, partyswitch,
                                                                             start_second_term, sameagency,
                                                                             z_pres_congdist, z_polarization,
                                                                             z_pres_agencydist, z_pres_fildist,
                                                                             divided, z_rollcalls, z_appts_total, central_staff_agency,
                                                                             prev_exec_pct, prev_ses_noncareer_pct_gap,
                                                                             looming_election, zmecompmedian))

  foo.newdata.unified.fil.total$z_pres_fildist <- dist.range[i]
  foo.newdata.unified.cong.total$z_pres_congdist <- dist.range[i]
  foo.newdata.divided.fil.total$z_pres_fildist <- dist.range[i]
  foo.newdata.divided.cong.total$z_pres_congdist <- dist.range[i]
  foo.newdata.unified.fil.total$divided <- 0
  foo.newdata.divided.fil.total$divided <- 1
  foo.newdata.unified.cong.total$divided <- 0
  foo.newdata.divided.cong.total$divided <- 1
  
  foo.ses.fil.preds.unified.total <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.unified.fil.total, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.unified.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.unified.fil.total))])
  
  foo.ses.cong.preds.unified.total <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.unified.cong.total, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.unified.cong.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.unified.cong.total))])
  
  foo.schedc.fil.preds.unified.total <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.unified.fil.total, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  
  foo.schedc.cong.preds.unified.total <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.unified.cong.total, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.cong.total))]
  
  foo.ses.fil.preds.divided.total <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.divided.fil.total, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.ses.cong.preds.divided.total <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.divided.cong.total, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.cong.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.cong.total))])
  
  foo.schedc.fil.preds.divided.total <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.divided.fil.total, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.divided.fil.total))]
  
  foo.schedc.cong.preds.divided.total <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.divided.cong.total, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.divided.cong.total))]
  
  ses.fil.preds.unified.total <- rbind(ses.fil.preds.unified.total, foo.ses.fil.preds.unified.total)
  ses.cong.preds.unified.total <- rbind(ses.cong.preds.unified.total, foo.ses.cong.preds.unified.total)
  schedc.fil.preds.unified.total <- rbind(schedc.fil.preds.unified.total, foo.schedc.fil.preds.unified.total)
  schedc.cong.preds.unified.total <- rbind(schedc.cong.preds.unified.total, foo.schedc.cong.preds.unified.total)
  ses.fil.preds.divided.total <- rbind(ses.fil.preds.divided.total, foo.ses.fil.preds.divided.total)
  ses.cong.preds.divided.total <- rbind(ses.cong.preds.divided.total, foo.ses.cong.preds.divided.total)
  schedc.fil.preds.divided.total <- rbind(schedc.fil.preds.divided.total, foo.schedc.fil.preds.divided.total)
  schedc.cong.preds.divided.total <- rbind(schedc.cong.preds.divided.total, foo.schedc.cong.preds.divided.total)
  
}

# partisan/ideology differences (-1 to 1 sd diff)
numappoint.diffs <- rbind(t(apply(diff(schedc.cong.preds.unified.total[c(26, 76),]),1,quantile, c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T)),
                          t(apply(diff(schedc.fil.preds.unified.total[c(26, 76),]),1,quantile, c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T)),
                          t(apply(diff(schedc.cong.preds.divided.total[c(26, 76),]),1,quantile, c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T)),
                          t(apply(diff(schedc.fil.preds.divided.total[c(26, 76),]),1,quantile, c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T)),
                          t(apply(diff(ses.cong.preds.unified.total[c(26, 76),]),1,quantile, c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T)),
                          t(apply(diff(ses.fil.preds.unified.total[c(26, 76),]),1,quantile, c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T)),
                          t(apply(diff(ses.cong.preds.divided.total[c(26, 76),]),1,quantile, c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T)),
                          t(apply(diff(ses.fil.preds.divided.total[c(26, 76),]),1,quantile, c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T)))

rownames(numappoint.diffs) <- c(1:8)
numappoint.diffs <- as.data.frame(numappoint.diffs)
colnames(numappoint.diffs) <- c("low95", "low90", "mid", "high90", "high95")
numappoint.diffs$divided <- c("Unified Government", "Unified Government",
                              "Divided Government", "Divided Government",
                              "Unified Government", "Unified Government",
                              "Divided Government", "Divided Government")

numappoint.diffs$pivot <- c("President-Senate Distance", "President-Filibuster Distance",
                            "President-Senate Distance", "President-Filibuster Distance",
                            "President-Senate Distance", "President-Filibuster Distance",
                            "President-Senate Distance", "President-Filibuster Distance")

numappoint.diffs$type <- c("Number of Schedule C Positions", "Number of Schedule C Positions", 
                           "Number of Schedule C Positions", "Number of Schedule C Positions",
                           "Number of Noncareer SES Executives","Number of Noncareer SES Executives",
                           "Number of Noncareer SES Executives", "Number of Noncareer SES Executives")

fig3 <- ggplot(numappoint.diffs, aes(x = divided, y = mid, shape = pivot, color = pivot)) + 
  geom_point(position = position_dodge(width = 0.25)) + 
  facet_wrap(~type, scales = "free_y", ncol = 1) +
  geom_errorbar(aes(ymin = low90, ymax = high90), width = 0, size = 0.75, position = position_dodge(width = 0.25)) + 
  geom_errorbar(aes(ymin = low95, ymax = high95), width = 0, position = position_dodge(width = 0.25)) + 
  geom_hline(yintercept = 0, linetype = "dotted") + 
  theme_bw() + 
  ylab('Estimated Change in Number of Indicated Position Type\n(Shifting Ideology from One SD Below Mean to One SD Above)\n') + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  xlab('\nPartisan Control of Government\n') + 
  scale_shape_discrete(name = "Ideological Measure") + 
  scale_color_grey(name = "Ideological Measure", start = 0.2, end = 0.6)

ggsave(fig3, file = "Figure-3.pdf", height = 6, width = 6)
ggsave(fig3, file = "Figure-3.png", height = 6, width = 6)



# now do this for the rest of the independent variables
ses.glmer.fil.preds <- ses.glmer.cong.preds <- schedc.glmer.fil.preds <- schedc.glmer.cong.preds <-
  data.frame(stdev = c("low", "high"),
             low95 = NA,
             low90 = NA,
             pred = NA,
             high90 = NA,
             high95 = NA)

  foo.newdata <- subset(apptypes_merged, agency != "All",
                        select = c(pres_switch, partyswitch, start_second_term, sameagency,
                                   z_pres_congdist, z_polarization, z_pres_agencydist, z_pres_fildist,
                                   divided, z_rollcalls, z_appts_total, central_staff_agency,
                                   prev_exec_pct, prev_ses_noncareer_pct_gap, looming_election, zmecompmedian))
  
  foo.newdata.pqexp.low <- foo.newdata.pqexp.high <- foo.newdata
  foo.newdata.pqexp.low$prev_exec_pct <-  mean(apptypes_merged$prev_exec_pct[as.numeric(rownames(foo.newdata.divided.fil.total))], na.rm = T) - 
    sd(apptypes_merged$prev_exec_pct[as.numeric(rownames(foo.newdata.divided.fil.total))], na.rm = T)
  foo.newdata.pqexp.high$prev_exec_pct <-  mean(apptypes_merged$prev_exec_pct[as.numeric(rownames(foo.newdata.divided.fil.total))], na.rm = T) + 
    sd(apptypes_merged$prev_exec_pct[as.numeric(rownames(foo.newdata.divided.fil.total))], na.rm = T)
  
  foo.ses.fil.preds.pqexp.low <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.pqexp.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.pqexp.low <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.pqexp.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.pqexp.low <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.pqexp.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.pqexp.low <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.pqexp.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  
  foo.ses.fil.preds.pqexp.high <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.pqexp.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.pqexp.high <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.pqexp.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.pqexp.high <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.pqexp.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.pqexp.high <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.pqexp.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  
  foo.newdata.noncareer.low <- foo.newdata.noncareer.high <- foo.newdata
  foo.newdata.noncareer.low$prev_ses_noncareer_pct_gap <-  mean(apptypes_merged$prev_ses_noncareer_pct_gap[as.numeric(rownames(foo.newdata.divided.fil.total))], na.rm = T) - 
    sd(apptypes_merged$prev_ses_noncareer_pct_gap[as.numeric(rownames(foo.newdata.divided.fil.total))], na.rm = T)
  foo.newdata.noncareer.high$prev_ses_noncareer_pct_gap <-  mean(apptypes_merged$prev_ses_noncareer_pct_gap[as.numeric(rownames(foo.newdata.divided.fil.total))], na.rm = T) + 
    sd(apptypes_merged$prev_ses_noncareer_pct_gap[as.numeric(rownames(foo.newdata.divided.fil.total))], na.rm = T)
  
  foo.ses.fil.preds.noncareer.low <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.noncareer.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.noncareer.low <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.noncareer.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  
  foo.ses.fil.preds.noncareer.high <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.noncareer.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.noncareer.high <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.noncareer.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.newdata.sameagency.low <- foo.newdata.sameagency.high <- foo.newdata
  foo.newdata.sameagency.low$sameagency <- -1
  foo.newdata.sameagency.high$sameagency <- 1
  
  foo.ses.fil.preds.sameagency.low <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.sameagency.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.sameagency.low <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.sameagency.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.sameagency.low <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.sameagency.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.sameagency.low <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.sameagency.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  
  foo.ses.fil.preds.sameagency.high <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.sameagency.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.sameagency.high <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.sameagency.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.sameagency.high <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.sameagency.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.sameagency.high <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.sameagency.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  
  foo.newdata.newpresident.low <- foo.newdata.newpresident.high <- foo.newdata
  foo.newdata.newpresident.low$pres_switch <- 0
  foo.newdata.newpresident.high$pres_switch <- 1
  
  foo.ses.fil.preds.newpresident.low <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.newpresident.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.newpresident.low <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.newpresident.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.newpresident.low <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.newpresident.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.newpresident.low <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.newpresident.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  
  foo.ses.fil.preds.newpresident.high <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.newpresident.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.newpresident.high <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.newpresident.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.newpresident.high <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.newpresident.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.newpresident.high <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.newpresident.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  
  foo.newdata.newparty.low <- foo.newdata.newparty.high <- foo.newdata
  foo.newdata.newparty.low$partyswitch <- 0
  foo.newdata.newparty.high$partyswitch <- 1
  
  foo.ses.fil.preds.newparty.low <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.newparty.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.newparty.low <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.newparty.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.newparty.low <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.newparty.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.newparty.low <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.newparty.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  
  foo.ses.fil.preds.newparty.high <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.newparty.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.newparty.high <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.newparty.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.newparty.high <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.newparty.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.newparty.high <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.newparty.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  
  foo.newdata.start_second_term.low <- foo.newdata.start_second_term.high <- foo.newdata
  foo.newdata.start_second_term.low$start_second_term <- 0
  foo.newdata.start_second_term.high$start_second_term <- 1
  
  foo.ses.fil.preds.start_second_term.low <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.start_second_term.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.start_second_term.low <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.start_second_term.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.start_second_term.low <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.start_second_term.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.start_second_term.low <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.start_second_term.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  
  foo.ses.fil.preds.start_second_term.high <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.start_second_term.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.start_second_term.high <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.start_second_term.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.start_second_term.high <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.start_second_term.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.start_second_term.high <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.start_second_term.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]

  foo.newdata.looming_election.low <- foo.newdata.looming_election.high <- foo.newdata
  foo.newdata.looming_election.low$looming_election <- 0
  foo.newdata.looming_election.high$looming_election <- 1
  
  foo.ses.fil.preds.looming_election.low <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.looming_election.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.looming_election.low <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.looming_election.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.looming_election.low <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.looming_election.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.looming_election.low <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.looming_election.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  
  foo.ses.fil.preds.looming_election.high <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.looming_election.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.looming_election.high <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.looming_election.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.looming_election.high <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.looming_election.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.looming_election.high <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.looming_election.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  
  foo.newdata.rollcalls.low <- foo.newdata.rollcalls.high <- foo.newdata
  foo.newdata.rollcalls.low$z_rollcalls <- -1
  foo.newdata.rollcalls.high$z_rollcalls <- 1
  
  foo.ses.fil.preds.rollcalls.low <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.rollcalls.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.rollcalls.low <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.rollcalls.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.rollcalls.low <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.rollcalls.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.rollcalls.low <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.rollcalls.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  
  foo.ses.fil.preds.rollcalls.high <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.rollcalls.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.rollcalls.high <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.rollcalls.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.rollcalls.high <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.rollcalls.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.rollcalls.high <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.rollcalls.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  
  
  foo.newdata.polarization.low <- foo.newdata.polarization.high <- foo.newdata
  foo.newdata.polarization.low$z_polarization <- -1
  foo.newdata.polarization.high$z_polarization <- 1
  
  foo.ses.fil.preds.polarization.low <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.polarization.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.polarization.low <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.polarization.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.polarization.low <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.polarization.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.polarization.low <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.polarization.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  
  foo.ses.fil.preds.polarization.high <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.polarization.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.polarization.high <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.polarization.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.polarization.high <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.polarization.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.polarization.high <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.polarization.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  
  
  foo.newdata.mecompmedian.low <- foo.newdata.mecompmedian.high <- foo.newdata
  foo.newdata.mecompmedian.low$zmecompmedian <- -1
  foo.newdata.mecompmedian.high$zmecompmedian <- 1
  
  foo.ses.fil.preds.mecompmedian.low <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.mecompmedian.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.mecompmedian.low <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.mecompmedian.low, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.mecompmedian.low <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.mecompmedian.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.mecompmedian.low <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.mecompmedian.low, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  
  foo.ses.fil.preds.mecompmedian.high <- plogis(predict(ses.glmer.fil.expertise, newdata = foo.newdata.mecompmedian.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  foo.ses.cong.preds.mecompmedian.high <- plogis(predict(ses.glmer.cong.expertise, newdata = foo.newdata.mecompmedian.high, re.form = NA))*
    (apptypes_merged$perm_ses_career[as.numeric(rownames(foo.newdata.divided.fil.total))] + 
       apptypes_merged$perm_ses_noncareer[as.numeric(rownames(foo.newdata.divided.fil.total))])
  
  foo.schedc.fil.preds.mecompmedian.high <- plogis(predict(schedc.glmer.fil.expertise, newdata = foo.newdata.mecompmedian.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]
  foo.schedc.cong.preds.mecompmedian.high <- plogis(predict(schedc.glmer.cong.expertise, newdata = foo.newdata.mecompmedian.high, re.form = NA))*
    apptypes_merged$total[as.numeric(rownames(foo.newdata.unified.fil.total))]  
  
    
control.mfxs <- data.frame(diff = rbind(quantile(foo.schedc.cong.preds.pqexp.high - foo.schedc.cong.preds.pqexp.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.cong.preds.pqexp.high - foo.ses.cong.preds.pqexp.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.schedc.fil.preds.pqexp.high - foo.schedc.fil.preds.pqexp.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.fil.preds.pqexp.high - foo.ses.fil.preds.pqexp.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.schedc.cong.preds.polarization.high - foo.schedc.cong.preds.polarization.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.cong.preds.polarization.high - foo.ses.cong.preds.polarization.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.schedc.fil.preds.polarization.high - foo.schedc.fil.preds.polarization.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.fil.preds.polarization.high - foo.ses.fil.preds.polarization.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.schedc.cong.preds.sameagency.high - foo.schedc.cong.preds.sameagency.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.cong.preds.sameagency.high - foo.ses.cong.preds.sameagency.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.schedc.fil.preds.sameagency.high - foo.schedc.fil.preds.sameagency.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.fil.preds.sameagency.high - foo.ses.fil.preds.sameagency.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.schedc.cong.preds.newpresident.high - foo.schedc.cong.preds.newpresident.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.cong.preds.newpresident.high - foo.ses.cong.preds.newpresident.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.schedc.fil.preds.newpresident.high - foo.schedc.fil.preds.newpresident.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.fil.preds.newpresident.high - foo.ses.fil.preds.newpresident.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.schedc.cong.preds.newparty.high - foo.schedc.cong.preds.newparty.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.cong.preds.newparty.high - foo.ses.cong.preds.newparty.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.schedc.fil.preds.newparty.high - foo.schedc.fil.preds.newparty.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.fil.preds.newparty.high - foo.ses.fil.preds.newparty.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.schedc.cong.preds.start_second_term.high - foo.schedc.cong.preds.start_second_term.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.cong.preds.start_second_term.high - foo.ses.cong.preds.start_second_term.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.schedc.fil.preds.start_second_term.high - foo.schedc.fil.preds.start_second_term.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.fil.preds.start_second_term.high - foo.ses.fil.preds.start_second_term.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.schedc.cong.preds.looming_election.high - foo.schedc.cong.preds.looming_election.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.cong.preds.looming_election.high - foo.ses.cong.preds.looming_election.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.schedc.fil.preds.looming_election.high - foo.schedc.fil.preds.looming_election.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.fil.preds.looming_election.high - foo.ses.fil.preds.looming_election.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.schedc.cong.preds.rollcalls.high - foo.schedc.cong.preds.rollcalls.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.cong.preds.rollcalls.high - foo.ses.cong.preds.rollcalls.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.schedc.fil.preds.rollcalls.high - foo.schedc.fil.preds.rollcalls.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.fil.preds.rollcalls.high - foo.ses.fil.preds.rollcalls.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.schedc.cong.preds.mecompmedian.high - foo.schedc.cong.preds.mecompmedian.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.cong.preds.mecompmedian.high - foo.ses.cong.preds.mecompmedian.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.schedc.fil.preds.mecompmedian.high - foo.schedc.fil.preds.mecompmedian.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.fil.preds.mecompmedian.high - foo.ses.fil.preds.mecompmedian.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.cong.preds.noncareer.high - foo.ses.cong.preds.noncareer.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T),
                          quantile(foo.ses.fil.preds.noncareer.high - foo.ses.fil.preds.noncareer.low, 
                                   c(0.025, 0.05, 0.5, 0.95, 0.975), na.rm = T)))
  
control.mfxs$variable <- c(rep("Previous Quarter Executive Schedule Percentage", 4),
                           rep("Congressional Polarization", 4),
                           rep("President-Agency Alignment", 4),
                           rep("New President", 4),
                           rep("New Presidential Party", 4),
                           rep("Presidential Reelection", 4),
                           rep("Upcoming Presidential Election", 4),
                           rep("Quarterly Senate Roll Calls", 4),
                           rep("Managerial Expertise", 4),
                           rep("Previous Quarter Noncareer SES Gap", 2))  
  
control.mfxs$variable <- factor(control.mfxs$variable, levels = c("Congressional Polarization", 
                                                                  "Managerial Expertise", 
                                                                  "President-Agency Alignment", 
                                                                  #"Central Staffing Agency",
                                                                  "Quarterly Senate Roll Calls",
                                                                  "Presidential Reelection", 
                                                                  "Upcoming Presidential Election", 
                                                                  "New President",
                                                                  "New Presidential Party", 
                                                                  "Previous Quarter Executive Schedule Percentage",
                                                                  "Previous Quarter Noncareer SES Gap"))

control.mfxs$pivot <- c(rep(c("President-Senate Distance", "President-Senate Distance",
                              "President-Filibuster Distance", "President-Filibuster Distance"), 9),
                        c("President-Senate Distance", "President-Filibuster Distance"))

control.mfxs$type <- c(rep(c("Number of Schedule C Positions", "Number of Noncareer SES Executives"), 18),
                        c("Number of Noncareer SES Executives", "Number of Noncareer SES Executives"))

colnames(control.mfxs)[1:5] <- c("low95", "low90", "mid", "high90", "high95")



fig5 <- ggplot(control.mfxs, aes(x = variable, y = mid, shape = pivot, color = pivot)) + 
  geom_point(position = position_dodge(width = 0.25)) + 
  facet_wrap(~type, scales = "free_y", ncol = 1) +
  geom_errorbar(aes(ymin = low90, ymax = high90), width = 0, size = 0.75, position = position_dodge(width = 0.25)) + 
  geom_errorbar(aes(ymin = low95, ymax = high95), width = 0, position = position_dodge(width = 0.25)) + 
  geom_hline(yintercept = 0, linetype = "dotted") + 
  theme_bw() + 
  ylab('Estimated Change in Number of Indicated Position Type\n') + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  xlab('\nIndependent Variable\n') + 
  scale_shape_discrete(name = "Ideological Measure") + 
  scale_y_continuous(trans = "asinh", breaks = c(-100, -10, -1, 0, 1, 10, 100, 1000)) + 
  scale_color_grey(name = "Ideological Measure", start = 0.2, end = 0.6)

ggsave(fig5, file = "Figure-5.pdf", height = 7, width = 10)
ggsave(fig5, file = "Figure-5.png", height = 7, width = 10)



